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ABSTRACT 

In a serie of three papers, the dynamical interplay between environments and dark mat- 
ter haloes is investigated, while focussing on the dynamical flows through the virtual virial 
sphere. It relies on both cosmological simulations, to constrains the environments, and an ex- 
tension to the classical m atrix method to derive the responses of the halo. A companion paper 
dPichon & Aubertl (E)06 ; ), paper I) showed how perturbation theory allows to propagate the 
statistical properties of the environment to an ensemble description of the dynamical response 
of the embedded halo. The current paper focuses on the statistical characterisation of the en- 
vironments surrounding haloes, using a set of large scale simulations; the large statistic of 
environments presented here allows to put quantitative and statistically significant constrains 
on the properties of flows accreted by haloes. 

The description chosen in this paper relies on a "fluid" halocentric representation. The 
interactions between the halo and its environment are investigated in terms of a time dependent 
external tidal field and a source term characterizing the infall. The former accounts for fly bys 
and interlopers. The latter stands for the distribution function of the matter accreted through 
the virial sphere. The method of separation of variables is used to decouple the temporal 
evolution of these two quantities from their angular and velocity dependence by means of 
projection on a 5D basis. 

It is shown how the flux densities of mass, momentum and energy can provide an al- 
ternative description to the 5D projection of the source. Such a description is well suited to 
re-generate synthetic time-lines of accretion which are consistent with environments found in 
simulations as discussed in appendix. The method leading to the measurements of these quan- 
tities in simulations is presented in details and applied to 15000 haloes, with masses between 
5 • 10 12 Mq and 10 14 M Q evolving between z = 1 and z = 0. The influence of resolution, 
class of mass, and selection biases are investigated with higher resolution simulations. The 
emphasis is put on the one and two-point statistics of the tidal field, and of the flux density of 
mass, while the full characterization of the other fields is postponed to paper III. 

The net accretion at the virial radius is found to decrease with time. This decline results 
from both an absolute decrease of infall and from a growing contribution of outflows. Infall 
is found to be mainly radial and occurring at velocities ~ 0.75 times the virial velocity. 
Outflows are also detected through the virial sphere and occur at lower velocities ~ 0.6V C on 
more circular orbits. The external tidal field is found to be strongly quadrupolar and mostly 
stationnary, possibly reflecting the distribution of matter in the halo's near environment. The 
coherence time of the small scale fluctuations of the potential hints a possible anisotropic 
distribution of accreted satellites. The flux density of mass on the virial sphere appears to 
be more clustered than the potential while the shape of its angular power spectrum seems 
stationnary. Most of these results are tabulated with simple fitting laws and are found to be 
consistent with published work which rely on a description of accretion in terms of satellites. 



Key words: Galaxies: formation, kinematics and dynamics. Cosmology. Methods: N-Body 
simulations. 
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1 GALAXIES IN THEIR ENVIRONMENT 

Examples of galaxies interacting with their environments are nu- 
merous. The antennae, the cartwheel galaxy, M51 are among the 
most famous ones. One of our closest neighbor, M31, exhibits a gi- 
ant stellar stream which m ay be associated with its satellites (e.g. 
iMcConnachie et all 120031) ). Even the Milky Way shows relics of 
past interactions wit h material coming from the outskirts, such as 
the Sagittarius dwarf llbata et al]ll995l) ). It appears clearly that the 
evolution of galactic systems cannot be understood only by consid- 
ering their internal properties but also by taking into account their 
environment. From a dynamical point of view it is still not clear 
for example if spirals in galaxies are induced by intrinsic un sta- 
ble modes fe.g. lLvnden-Bell & Kalnai si i 19721) . iKalnai sll 1 9771) ') or 
if they are du e to gravitational interacti ons with satellites or other 
galaxies (e.g. iToomre & Toomrel ll97 2 |)). Similarly, normal mod e 
theories of warps have be en proposed Isparke & Casertand Jl988l) . 
iHunter & Toomrel Jl969h ) but failed to reproduce l ong-lived warps 
in a live halo for example (e.g. lBinnev et a lj (1998])). Since warped 
galaxi es are likely to have companions iReshetnikov & Combesl 
( 1998)), it is nat ural to suggest sat e llite tidal forcing as a generating 
mechanism (e.g. lWeinber d l 1 998h . lTsu"chiviJ <2002l) ). Another pos- 
sibil ity is angular momentum misalignment of inf aili ng material 
(e.g. lOstriker & Binnevl Jl989h . ljiang & Binnevl <1999l) ). The exis- 
tenc e of the thick disk may also b e explained by past s mall mergers 
(e.g.lVelazquez & White! Il999h . l0uinn et alj fl993l) . I Walker etaH 
( 1996)). Conversely, very thin disk put serious constraints on the 
amplitude of the interactions they may have experienced in the past. 

On a larger scale, dark matter halos are built in a hierar- 
chical fashion within the CDM model. Some of the most seri- 
ous challenges that these models are now facing (t he over pro- 
ductio n of dwarf gal a xies i n the local group (e.g. iMoore et alj 
Jl999h . iKIypin et alj Jl999h) . the c u spide crisis o f NFW-like 
halos (e.g. iFlores & Primackl Jl994l) . iMoord <1994l) ). the over- 
cooling problem and the mo mentum crisis for galactic disks (e.g. 
iNavarro & Steinmetd 1 19971) ) occur at these scales; it is therefore 
important to study the effects of the cosmological paradigm on the 
evolution of galaxies in order to address these issues. 

In fact, the properties of galaxies naturally presen t correlations 
with their environments. For example iTormed 119971) showed that 
the shape of halos tends to be aligned with the surrounding satel- 
lites distribution. Also, h alos spin is sensitive to recently accrete d 
angul ar momentum (e.g. Ivan Haarlem & van deWevgaert|jl993l) . 
lAubert et alj 120041) ). More generally, halos inherit the properties 
of their progenitors. 

At this point, a question naturally arises, namely 'what is the 
dynamical response of a galactic system (halo+disk) to its environ- 
ment ?'. One way to address this issue is to compute hi gh resolution 
simula t ions of galax i es into a given environm ent (e.g. lAbadi et alj 
l2003l) . lKnebe et alj <2004 . lGill et alj Exil ). However, if one is 
interested in reproducing the variety of dynamical responses of 
galaxies to various environments, the use of such simulations be- 
comes rapidly tedious. An alternative way to investigate this topic 
is presented here, which should complement both high resolution 
simulations and large cosmological simulations. In a serie of three 
papers, an hybrid approach is presented to investigate the interplay 
between environments and haloes. It relies on both cosmological 
simulations (to constrains the environments) and on a straightfor- 
ward extension of the classical tools of galactic dyna mics (to de- 
rive th e haloes response). A companion paper IPichon & Aubertl 
fcOOd) . Paper I hereafter, in press) describes the analytic theory 
which allow to assess the dynamics of haloes in the open, secular 



and non-linear regimes. The purpose of the current paper is to set 
out a framework in which to describe statistically the environments 
of haloes and pre sent results on the tidal fie ld and the flux density of 
matter. Paper III ( Aubert & Pichod 12006I) . in prep.) will conclude 
the complete description of the environments of dark haloes. 



1.1 Galactic infall as a cosmic boundary 

Clearly a number of problems concerning galactic evolution can 
only be tackled properly via a detailed statistical investigation. Let 
us briefly make an analogy to the cosmological growth of density 
fluctuations. Under certain assumptions, one can solve the equa- 
tions of_erohitionof tho se over-densities in an ex panding universe 
(e.g. lPeeblesI ll98(l) . see lBernardeau et all 120021) for an extensive 
review). Their statistical evolution due to gravitational clustering 
follows, given the statistical properties of the initial density field. 
For example, the power spectrum, P(z, k), may be computed for 
various primordial power spectra, P pr im(fc), and for various cos- 
mologies. In other words, the statistical properties of the initial 
conditions are propagated to a given redshift through an operator 
N given by the non linear dynamical equations of the clustering : 

P{ Z ,k) = K(P pI im{k),z). (1) 

In a similar way, how would the statistical properties of en- 
vironments be propagated to the dynamical properties of galactic 
systems? This is clearly a daunting task : the previous analogy with 
the cosmological growth of perturbation is restricted to its princi- 
ple. For example the assumption of a uniform and cold initial state 
cannot be sustained for galaxies and halos. While spatial isotropy 
is clearly not satisfied by disks, and hot, possibly triaxial halos, the 
velocity tensor of galaxies may also be anisotropic. Environments 
also share these inhomogeneous and anisotropic features since they 
are also the product of gravitational clustering and cannot be sim- 
ply described as Gaussian fields. These boundary conditions are 
not pure 'initial conditions' since they evolve with time and in a 
non stationary manner (e.g. the accretion rate decreases with time). 
A whole range of mass must be taken into account, each with differ- 
ent statistical properties. Finally, trajectories cannot be considered 
as ballistic (even in the linear regime) and must be integrated over 
long periods. Notwithstanding the above specificities of the galactic 
framework, two questions have to be answered: 

(i) What is the 'galactic' equivalent of P pl -i m (fc), i.e. how to de- 
scribe statistically the boundary conditions ? 

(ii) What is the 'galactic' equivalent of i.e. how to describe 
the inner galactic dynamics ? 

The second point is discussed extensively in lPichon & Aubertl 
120061) and is briefly summarized in section [2] In that paper, it 
is shown how a perturbative theory can describe the dynamics of 
halo es which experience both accretion and tidal interactions (see 
also lAubert etall <2004 ). Within this formalism, the environment 
is described by the external gravitational potential and a source 
function. The former describes fly-bys and the tidal field of neigh- 
bouring large scale structures. The latter describes the flows of 
dark matter, i.e. the exchanges of material between the halo and 
the 'inter-halo' medium. The knowledge of these two quantities 
fully characterizes the boundary condition. The focus here is on 
well formed halos which do not undergo major merger between 
2 = 1 and z — 0. This bias is consistent with a galactocentric 
description in which a perturbative description of the inner dynam- 
ics is appropriate and equal mass mergers are explicitly ignored. 
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As briefly explained in section[2] this formalism provides a link be- 
tween the statistical properties of environments to the statistical dis- 
tributions of the responses of haloes: this link is referred to as sta- 
tistical propagation. In this manner, the distribution of haloes dy- 
namical state can be directly inferred from the statistical properties 
of environments, without relying on the follow up of individual in- 
teracting haloes. The observed distributions of dynamical features 
provides informations on the cosmic boundaries which influence 
haloes. This, together w ith the perturbative formalism described in 
IPichon & Aubertl l2006t) should allow us to address statistically the 
recurrent 'nurture or nature' problem of structure formation within 
galactic systems. 

This statistical formalism is complementary to methods based 
on merger trees (which also coupl e environment and inner prop - 
erties of galact i c systems see e.g. iKauffmann & White! Jl993J> . 
iRoukema et"ail <1997l)ISomerville & Kolattl <1999alV ). These 'ana- 
lytic' or 'semi-analytic' models, with presciption for the baryons 
contained in haloes, angular momentum transfer, cooling and star 
formation may pr edict properties o f galaxies given their forma- 
tion history (e.g. ICole et all 11994) ). This history may be pro- 
vid ed analytica l ly using extended Press-Sch echter formalism (see 
e.g.|Bondrt_alJ (l9^lb,|Lacev&Col3 jl993|) ) or us ing simulations 
(e.g. IKauffmann et alj7l999MBenson et all EoOll) ). Even tough 
this technic now extends its field of application to subhaloes (see 
e.g. iBlaizot et alj 120061) ). it remains somewhat limited for the 
purpose dynamical applications. These require a detailed descrip- 
tion of the geometrical configuration of the perturbations, and of 
the dynamical response of the halo. Both of these are difficult 
to reduce to simple recipes. Conversely, full analytic theories of 
the inner dynamics of interac ti ng haloes were d e velopped in e.g . 
iTremaine & Weinberg! <1984l) . IWeinberd <1998l) . iMuralil <19991) . 
Relying on the matrix method, these theories do take properly into 
account the resonant processes that occur when the halo is per- 
turbed by an external potential. However, they usually do not ac- 
count for the perturbations induced by the accretion of matter, while 
these authors generally considered test-cases where a halo respo nds 
to a given configuration (or statistics, see e.g. IWeinberdl200lh of 
perturbations. Paper I extended these theories to open stellar sys- 
tems and , while relying on numerical simulations to constrain the 
environments, it reformulated them in terms of the statistics of the 
inner dynamics of a representative population of haloes. 

Paper I presents a list of possible applications. For instance, 
gravitational lensing by halos is affected by inner density fluctu- 
ations, which are induced by the halo's environment : hence the 
statistics of the lensing signal is be related to the statistics of halo's 
perturbations, therefore to the cosmological growth of structures. 
Paper I showed how this approach could be extended to other ob- 
servables, such as X-Ray temperature maps, SZ surveys or direct 
detection of dark matter. Statistical propagation allows to relate 
cosmology to the inner properties of cluster and galaxies. Con- 
versely, it should be possible to show if the perturbations measured 
in simulations are consistent with a secular drift toward a universal 
profile of haloes. Closer to us, the correlation of the numerous arte- 
facts of past accretion in the local group, such as streams or tidal 
tails, can be understood in terms of statistics of environments. All 
processes which depend critically on the geometry of the interac- 
tions may be tackled in this framework 1 . 



1 However, all departure to angular isot ropy on the sphere w ill be ignored 
here (in contrast to what was stressed in lAubert et alj |2004|) ), and its im- 
plications will be postponed to the discussions in sectionlS! 



The statistical propagation relies on the knowledge of the 
properties of the environment and is stated by the point (i) men- 
tioned above. This question is investigated the current paper by us- 
ing a large set of simulations, where each halo provides a realiza- 
tion of the environment. From this large ensemble of interacting 
haloes, the aim is to extract the global properties of their "cosmic 
neighborhood". Such a task requires an appropriate description of 
the source and the surrounding tidal field. It is the purpose of this 
work to implement this description which should both provide in- 
sights on the generic properties of cosmic environments and be use- 
ful in a 'dynamical' context. Specifically, a method is presented to 
constrain the exchanges between the halo and its neighborhood, 
via the properties of accretion and potential measured on the virial 
sphere. The advantages, specificities and caveats (and the methods 
implemented to overcome them) provided by this halocentric ap- 
proach will be presented in this paper. 

As shown in the following sections, the source function is 
given by the phase space distribution function (DF hereafter) of the 
adverted material. As a consequence, its full characterization is a 
complex task since it involves sampling a 5 dimensional space and 
relies on the projection of its DF on a suitable 5D basis. In particu- 
lar, it is shown how such a description can be used to constrain the 
kinematic properties of accretion by dark matter haloes in cosmo- 
logical simulations. The detailed statistical characterization of the 
higher moments of the source is postponed to a paper III. An alter- 
native description of the source is also presented; it relies on flux 
densities through the virial sphere, i.e. the moments of the source 
DF. Even though it is less suited to the dynamical propagation, this 
alternative description is easier to achieve numerically and to in- 
terpret physically. In particular, it illustrates how the source term 
may be characterized statistically via its moments. The link be- 
tween these flux densities and the 5D projection of the source is 
discussed together with the one and two-point statistics of the flux 
densities of mass through haloes in simulations. Also, the mean of 
reprojecting the effect of the external gravitational potential inside 
the halo (through Gauss's theorem) while knowing its properties on 
the virial sphere are discussed. The potential's one and two points 
statistics are also investigated around simulated haloes and inter- 
preted. 

Finally (Section|FJ provide means of regenerating such flows 
ab initio from its tabulated statistical properties. Such tools yields a 
way to embed idealized simulations of galaxies into realistics cos- 
mological environments. 

The outline of the paper is the following: Section |2| presents 
briefly the dynamics of open collisionless systems and states the 
principle of statistical propagation. Section [3] presents the proce- 
dure used to compute the source term, and illustrates its imple- 
mentation on a given halo. The simulations and the corresponding 
selection biases of our sample are then described in section l4~2l 
Section [5] and |6| present the statistical measurements for one and 
two point statistics respectively. SectionQdraws a global picture of 
galactic infall on L* galaxies, while a discussion and conclusions 
follow in section|8| 

Among the different results described in this paper, the reader 
will find : 

• a statistical description of the external gravitational field felt 
by halos : the potential is found to be quadrupolar and stationnary. 

• a study of the evolution of accretion : accrection by dark mat- 
ter haloes decrease with time, while the outflows becomes more 
significant at recent times. 

• constrains on the trajectory of infalling material : accretion is 
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found to be essentially radial, while outflows are found to be more 
circular. 

• results on the two-point statistics of the external potential 
measured on the viral sphere : the potential provide hints of an 
anisotropic perturbation of the halo. 

• results on the two-point statistic of the accretion's distribution 
on the virial sphere: accretion is dominated by small scales fluctu- 
ations and has a shorter coherence than the external gravitational 
field. 



2 DYNAMICS OF OPEN COLLISION-LESS SYSTEMS 

The exchanges occurring between a halo and its environment can 
be characterized in several ways. One of the classical method in- 
volves building a merger tree where the whole history of forma- 
tion of a halo is expressed in terms of global properties o f its pro - 
genitors (e.g.lLacev & Cole] Jl993l) . lKauffmann & White] <1993bl) . 
ISomerville & KolattHl999bl) ). While well suited to study the evolu- 
tion of those characteristics, it cannot be directly applied to predict 
in details the halos' inner dynamic because of the lack of spatial in- 
formation on these interactions. One could track the whole (six di- 
mensional) phase-space history of all the progenitors, but not only 
would it be difficult to store in practice it would also not give in- 
formation on the influence of large scale structures through their 
gravita tional potential. In the present paper, following lAubert et alj 
120041) . it is suggested to measure the relevant quantities on a sur- 
face at the interface between the halo and the intergalactic medium. 
Accretion is described as a flux of particles through the halos' ex- 
ternal boundaries. 

This se ction presents an extension o f the formali s m dev el- 
oped by e.g. iTremaine & Weinberg! 1 19841) and iMuralil 1 19991) to 
open spherical collisionless systems. The dynamics of a dark matter 
spherical halo is obtained by solving the collisionless Boltzmann 
equation coupled with the Poisson equation : 



dtF + v • d r F - W • d v F = 0, 



A* = 4ttG / d 3 vF(v) 



(2) 
(3) 



where F(r, v, t) is the system's distribution function coupled to 
ty(r, t) = ip + ip e , the total gravitational potential (self-gravitating 
+ external perturbation). Note that, in a somewhat unconventional 
manner, tfj e refers here to the external potential, i.e. the tidal poten- 
tial created by the perturbations outside the boundary. The gravi- 
tational field of incoming particles is accounted for by the source 
term. Equation J2j coupled with Hamilton's equations is a conser- 
vation equation: 



d t F + V(uf) = 0, 



(4) 



where u = (v,-V$) and V = (<9 r , d v ). As a consequence, con- 
sidering a 'source of material' described by / e (w) located on a 
surface S(w) implies: 

8 t F + V(uF) = -&(S(w) - o)u ■ T^§r/e(w), (5) 

where w = (r, v) describes the phase space and a defines the sur- 
face boundary of the studied system (here So stands for the Dirac 
delta function). If this boundary is defined as a spherical s urface 
with radius R, then equation (f5j becomes (e.g. lAppelH2002l) ) : 



d t F + V(u.F) 



-8u(r - i?)l> r / e (w) 
-5 D (r-R)s e (w). 



(6) 
(7) 



The function s e will be hereafter referred to as the 'source' func- 
tion. Formally the r.h.s. of equation can be seen as an additional 
local rate of change of the system's distribution function. Note that 
equation Q involves the external potential, t/> e , via u. 



2.1 Moments of the source term 

Integrating equation over velocities leads to the mass conserva- 
tion relation: 

dtp + V(vp) = -Sn{r - R){pv r ) e = -<5 D (r - R)zu p , (8) 

where the source appears as an external flux density of matter 
(pTv)e or vj p . Taking the next moment of equation Q leads to 
the Euler- Jeans Equation: 



dtpv + V • (pvv) + pV* = -fo(r - R)(pv r v) e 



(9) 



where the source adds a flux density of momentum, cc pv , to the 
conventional Jeans equation. Taking the successive moments of 
equation Q W1 U genetically include a new term in the resulting 
equations. 



2.2 Propagating the dynamics 

Following ITremaine & Weinberg! Il984t) . equation can be 
solved along with Poisson's equation in the regime of small per- 
turbations. In the spirit of paper let us define the system's environ- 
ment by the external perturbative potential ip e (r, t), and the source 
s B (r, v, t). Given this environment, the system's linear potential re- 
sponse i/j(r, t) can be computed. Writing the following expansions: 



V(r,t) = ^a n (^ [nl (r), 

n 

r(r,t) = ^& n (^ [nl (r), 

n 

s e (r,v,t) = ^c n (t)0 lnl (r,v), 



(10) 
(11) 
(12) 



where ^ [n, (r,v) and V [n] ( r ) 

are suitable basis functions, and 
solving the Boltzma nn and Poisson's equation for a n , one finds 
lAubertetal]j2004 ) : 

/OO 
drK(r - t) ■ [a(r) + b(r)] + H(r - t) ■ c(r). (13) 
- OO 

The kernels K and H are functions of the equilibrium state distri- 
bution function, Fo, and of the two basis, 0' n '(r, v), and (r) 
only (see paper I). As a consequence, they may be computed once 
and for all for a given equilibrium model. Since the basis function, 
V>'"', can be customized to the NFW-like profile of dark matter ha- 
los, it solves consistently and efficiently the coupled dynamical and 
field equations so long as the entering fluxes of dark matter amounts 
to a small perturbation in mass compared to the underlying equilib- 
rium. 

Assuming linearity and knowledge of K and H, one can see 
that the properties of the environments (through b and c) are prop- 
agated exactly to the inner dynamical properties of collisionless 
systems. Note in particular that the whole phase space response of 
the halo follow from the knowledge of a. For example, taking the 
temporal Fourier transform of equation 1 1 3i . the cross correlation 
matrix is easily deduced: 



> 



((I-K) • K-b + H- 
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[K • b + H ■ c] * • (I — K 



(14) 



where x = x(w) is the Fourier Transform of x(t). The environ- 
ment's two-points statistic, via (b • b* ), (c-c* ) and (b • c* ), 
modifies the correlation of the response of the inner halo. 

Linear dynamics do not take in account the effects on the per- 
turbation induced by dynamical friction. More generally the damp- 
ing of incoming fluxes will ultimately require non linear dynamics 
(since the relative temporal phases of the infall do matter in that 
context). It is also assumed in equation J 1 3 1 that the incoming mate- 
rial does not modify the equilibrium state of the system. The secular 
evolution of the system should also be ultimately taken in account, 
through a quasi-linear theory for example (see e.g. paper I). 

Let us emphasize that, since the addressed problem is lin- 
ear, the response, equation 1 1 3i . can be recast into a formulation 
which only involve an 'external potential', namely the sum of %j) e 
and the potential created by the entering particles described by s E . 
While formally simpler at the linear deterministic level, this alter- 
native formulation does not translate well non linearly or statisti- 
cally (since it would require the full knowledge of the perturbation 
everywhere in space in a manner which is dependent upon the inner 
structure of the halo). 

In the following sections, our aim is to describe how the two 
fields i/ ,e ( r i t) an d s e (r, v, t) can be extracted from haloes in cos- 
mological simulations. Then it will be shown how to characterize 
their statistical properties as a function of time via their expansion 
coefficients , b n (t) and c n (t). 

2.3 Convention and notations 

In what follows, let us characterize the properties of 2 fields, either 
angularly, kinematically, statistically or temporally, or any combi- 
nation (for various classes of masses). For a given field, X, let us 
introduce the following notations for clarity : 



IeI I X(0, (/>)dsm(e)d(j> , 



(15) 



which represents the angular average of X over the sphere. Alter- 
natively, let us define the temporal average over At as: 



-= AT I X(t)AU 



Finally let us define the ensemble average as 



(X) = / XF(X)dX = E{X}. 



(16) 



(17) 



where T is the density probability distribution of X. Here E{X} 
stands for the expectation of X. In practice, in section l4~2l an esti- 
mator for ensemble average of X measured for N halos is given : 2 



N 



(18) 



The underlying probability distributions are sometimes very 
skewed (when, e.g. corresponding to strong or weak accretion event 
around massive or smaller halos), which requires special care when 

2 An alternative would be to weight the sum by the relative number of 
objects in each halo, hereby down weighting light halos. It is found that this 
alternative estimator did not significantly affect our measurements 



attempting to define statistical trends. Hence let us also define < 
as the mode (or most probable value) of the fitting distribution of 
T. 

All external quantities, (flux densities, potential etc) will gen- 
erally be labeled as X e . Let us introduce moments of the source 
over velocities, which correspond to flux densities, noted vox and 
their corresponding fluxes, noted $x. Table 13.31 gives a list of 
such flux densities, flux pairs. Finally the harmonic transform of 
the field, X, will be written as af m and its corresponding power 
spectrum , while the parameters relative to fitting the statistics 
of the field will be written as qx ■ Note that the contrast of the field, 
X was also introduced as 



X -X 

(X) 



(19) 



and its corresponding harmonic transform, a em . A summary of all 
the notations can be found in Appendix IHI 



3 THE SOURCE OF INFALL 

Let us first describe our strategy to fully characterize the source 
of cosmic infall at the virial radius via collisionless dark matter 
simulations, and enumerate the corresponding biases. In particular 
let us illustrate our procedure on a template halo. 



3.1 Describing the source 

As argued in section|2| computing the response of an open system 
to infalling material requires the knowledge of the source function, 
s e (r, v,i). Given the particles accreted by a halo, one possibil- 
ity involves storing those phase-space properties for all particles. 
While feasible for a limited number of halos, this task would be- 
come rapidly intractable for our large number of simulations. In or- 
der to compress the information, the accreted distribution function 
is projected here on a basis of function, following equation J 121 . 

Since the measurement is carried at a fixed radius, the phase 
space is reduced from six to five degrees of freedom: two for the an- 
gular position on the sphere, described by two angles, (9, (f>) = fl, 
and three for the velocity space described in spherical coordinates 
by (v, Ti, T2) = (v, r), where v is the velocity modulus and T are 
the two angles describing its orientation (see figure^. The angle 
Ti indicates how radial is the velocity, with Ti > n/2 for infalling 
dark matter and Ti < n/2 for outflows. T2 indicates the orienta- 
tion of the infall's tangential motion. 

Recall that the two fields, c n (hence s e (tt, T, v, t)) and fe n 
(hence ip e (Cl, t)) are respectively five dimensional and two dimen- 
sional (as a function of mass and time). Note also that both s e and 
ip e are statistically stationary with respect to ft, while s e is par- 
tially isotropic and not stationary with respect to T; neither ip e and 
s e are stationary with respect to cosmic time. 

3.1.1 Harmonic expansion of the incoming fluxes 

The fi and T dependence are naturally projected on a basis of 
spherical harmonics, Yi im (Sl) and Y e i im i(T). The velocity am- 
plitude dependence is projected on a basis of Gaussian functions, 
g a (v) with means /j, a and a given r.m.s. a. One can write: 



' [nl (r,v) = Y m (Sl)Y m ,{T)g a (v), 



(20) 



where n = (I, m, a, £' , m') = (m, a, m'). The expansion coeffi- 
cients, c^Jf (t) are given by : 
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Figure 1. The angles f2 and V. The dot indicates the position of the particle 
on the sphere. The dashed ellipse represents the plane which contains both 
the particle and the sphere center. X and Z are arbitrary directions defined 
by the simulation's box. (8, (ft) = CI are the particle's angular coordinates 
on the sphere. V = (r± , 1^2 ) define the orientation of the particle's velocity 
vector (shown as an arrow). 

Cm'enC*) = • S™/) a . (21) 

where: 

(sSO/9 = / dmrdvv 2 g (v)Y^(fl)Y m ,(fl)s e (v,fl,T,t),(22) 
given 

Ga,g = J dvv 2 g a (v)gp(v). (23) 

Note that the expansion defined in equation 1 1 21 where the coef- 
ficients are given by equation 12 1 1 involves 5 subscripts spanning 
the 5 dimensional phase space, while the expansion in equation i 101 
only involve 3 subscripts. This description of the source term is re- 
duced to a set of coefficients which depends on time only. Further- 
more this procedure requires to parse the particles only once, and 
all the momenta (e.g. mass flux density, PDF of impact parameter 
...) of the source terms can be computed directly from these coeffi- 
cients. As a consequence, the statistics of momenta follow linearly 
from the statistics of coefficients only, as shown in section l4~2l 

3.1.2 Harmonic expansion of the external potential 

Let us call b' im (t) the harmonic coefficients of t he expansion of the 
external potential on the virial sphere. Following lMuraliHl999T) .let 
us expand the potential over the biorthogonal basis, (ttf, m , <4i' m ), 
so that 

= ^6nW W (r), (24) 

n 

where V [n] ( r ) = Y e m {fl)u^ m (r) . The first equality in equa- 
tion 1241 corresponds to the inner solution of the three dimensional 



potential whose boundary condition is given by Y^(fl)b' lm on the 
sphere of radius ifeoo (defined below). Since the basis is biorthog- 
onal, it follows that 




It is therefore straightforward to recover the coefficient of the 3D 
external potential from that of the potential on the sphere. 



3.2 From simulations to expansion coefficients 

Once a halo is detected, its outer 'boundary' is defined as a sphere 
centered on its center of mass with a radius, i?200 ( or Virial radius), 
defined implicitly by 3A//(47ri?2oo) = 200po. This choice of ra- 
dius is the result of a compromise between being a large distance 
to the halo center, to limit the contribution of halo's inner material 
to fluxes, and being still close enough to the halo's border, to limit 
the simulation's fraction to be processed and avoid contributions 
of fly by objects. Let us emphasize that several definition of the 
virial radius can be found in the litterature, involving e.g. the crit- 
ical density, or a different contrast factor, where the latter may or 
may not depend on the cosmology. Hence, one should keep in mind 
that all the quantitative results presented in this article depend on 
our specific choice of a definition. 

The time evolution of accretion is measured backwards in time 
by following the biggest progenitor of each halo detected at redshift 
2 = 0. The positions and velocities of particles passing through the 
virial sphere between snapshots are then stored. All positions are 
measured relative to the biggest progenitor center of mass while 
velocities are measured relative to its average velocity, for each red- 
shift z. In one of the simulation described below, the total comoving 
drift distance of the center of mass was compared to the distance 
between the halo's position between z=l and z=0. The haloes were 
chosen to satisfy the criteria described in 14.21 It is found that the 
scattering of the motion of the center of mass represents less than 
10 % of the distance covered in 8 Gyrs. The center of mass of the 
biggest progenitor seems stable enough to be a reference. 

Sticking to the previous definition of -R200 would imply a 
changing outer boundary and an 'inertial' flux through a moving 
surface would have to be taken in account. To overcome this ef- 
fect, the sphere was kept constant in time at a radius equals to 
-R200 (z — 0). This choice corresponds to a reasonable approxima- 
tion since the actual virial radius does not change significantly with 
time between z < 2 for a reasonably smooth accretion history. As 
shown in figure|2|the virial radius at z = 1 is only 20% smaller than 
-R200 measured at 2 = 0. Larger haloes have larger variations but 
the median value of the difference between the two radii remains 
smaller than 30% for final masses smaller than 1O 14 M0. Finally, 
measurements were done using physical coordinates (and not co- 
moving coordinates). These choices were partly guided by the fact 
that they simplify future applications of these results to the inner 
dynamic of the halos (see paper I). 



3.2.1 Sampling on the sphere 

As shown in section|2| the source function s e reads 

s e = /(r, v,t)v r = 6'v(r - n(t))5v(v - Vi(f))iv,<- (26) 
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Figure 2. Top: the distribution of the ratio between the virial radius mea- 
sured at z = 1 and at z = 0, Bottom: the distribution of iJ20o( z = 
l)/i?200 ( z = 0) as a function of the halo's final mass. Each point represent 
one halo. Symbols stand for the median value of i?20o(z = l)/R20o(z = 
0) in 6 different classes of masses. Bars stand for the interquartile. The two 
measurements were performed on 9023 haloes which satisfy the selection 
criteria defined in section l4!2l 



Switching to spherical coordinates leads to: 



E 



<5d(-R2oo - r l {t)) 5o(v - Vi(t)) 



r>2 

-"•200 



sinfii sin(Ti) 
where i is the particle index. Now: 



Vr,i{t), 



(27) 



(28) 



Vr,i5-D(R200 — Ti(t)) = 



k 

E 

k 



So(t — t200,k,i) 



Wk,iSl){t — t 2 00,k,i) , 



(29) 



where t 2 oo,k.i corresponds to the fc-th passage of the i-th particle 
through the virtual boundary -R200 (and tv,fc,» the corresponding 
radial velocity). In our conventions, the weight function Wk,i takes 
the value 1 if the particle is entering and - 1 if its exiting the virial 
sphere. Given that our time resolution is finite, let us consider a 
time interval AT around t and define the (temporal) average phase 
space flux density over AT: 
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Figure 3. The distribution of interpolated radial velocities v r of particles 
passing through the virial radius. Those particles were taken from the whole 
history of accretion of a typical halo (-R200 = 860 kpc, M(z = 0) = 
3 ■ lO ia M0). Entering particles (solid line) v r < while exiting particles 
(dashed line) have v r > 0, as it should be. 



1 pt+AT 

- {t) = / rfrs£(r) - (30) 

Equation 1281 becomes: 

N 

1 W = > o act. d2 r ~rF~\ — Wi ' k - ( 31 ) 

^ v 2 AT ■ iijoo sinfli sm(ri) 

i,k 

The simulations were sampled in time regularly in ln(z) (i.e. 
Aln(z) = constant). From z = 2 to z — 0.1, 23 snapshots 
were taken (and a z — snapshot was added to the sample). If At 
is small, the sum over k should mostly involve one passage, i.e. : 



AT 



1 5(v - 



s{v - v,) sjn-n,) c5(r-rj) 

sinfii sin(ri) 



Now these measurements only give access to (v, CI, T) at fixed red- 
shift, z, every varying Az. Consequently, these values need to be 
interpolated at the sought £200,1 approximated by : 

t 2 oo, ! = U(z n ) + t{ . Zn+l) ~ t{z ;\ (R2 00 - n{z n )). (33) 
ri(z n +i) — ri(z n ) 

Given these 'crossing' instant, the positions, r, and velocities, v, 
are also linearly interpolated. For instance, one gets for the x com- 
ponent of the velocity : 



i(feoo) = V Xt i(z„) + - 



,i(Zn + 1) 



Such an interpolation is not strictly self-consistent since a ballistic 
motion requires a constant velocity along the trajectory. The worst- 
case scenario would correspond to particles which have entered the 
virial sphere with an outflowing velocity vector and vice- versa. As 
a simple but important check, the distribution of interpolated radial 
velocities was plotted (see Fig. {5J). Those were computed from 
the whole history of accretion of a typical halo (-R200 = 860 kpc, 
M z= o = 3 • 1O 13 M0). The two types of particles (entering/exiting) 
are confined in their radial velocity plane: entering (resp. exiting) 
particles have negative (resp. positive) radial velocities. Velocities 
are correctly interpolated. It also means that our time steps are 
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Figure 4. The impact of time averaging on the measured scales of dark mat- 
ter passing through the sphere. On the left, time-position diagram of dark 
matter (black ellipses) as they pass through the sphere. Time integration is 
performed during AT (the two horizontal lines). On the right the accreted 
dark matter as seen on the sphere. A longer integration time increases the 
length scale of the incoming blob. If AT gets very large, different blobs 
may be seen as one (upper diagram). 



small enough to ensure a small variation of positions/velocities of 
particles, validating a posteriori our assumptions. A fraction of ex- 
iting particles do have a negative radial velocity but represent less 
than a few percent of the total population. For safety, those particles 
are rejected from the following analysis. 

One should note that the measured angular scales are sensi- 
tive to the time sampling (see Fig. J4j). Increasing the sampling 
time tends to increase the apparent size of objects as measured on 
the sphere. Since this increase depend on the shape or the orienta- 
tion of the objects, this effect cannot be simply time-averaged. As 
a consequence, a varying time step would induce a variation of typ- 
ical spatial scale. The interpolation given by equation I33i allows 
also for a constant time step resampling of the source term s_ e . All 
reference to the time average shall be dropped from now on. 

Given equation <3 It . computing the expansions coefficients of 
s e is straightforward : 



= G 



N 
-1 \ "* 



AT 



£?2 
"200 



(35) 



It is expected that t he above procedure is more accurate than the 
strategy presented in I Aubert et alj 120041) . where the flux densities 
were smoothed over a shell of finite thickness (-R200/IO). 

The harmonic expansion b' m of the external potential 
V> e (-R200, SI) is computed direclty from the positions of external 
particles (e.g. lMurali & Tremaind fl998l) ) : 



b'e,m(t) = 



i^V>/ m (fi j(i ))-52» 



21 + 1 ^ 

3 



/+!(*) 



(36) 



position of the j-th external particles. The quantities rj(t) and 
Slj(t) at time t are obtained by linear interpolation between two 
snapshots. Using equation 1241 - if} e (r < J?200,^) can be recon- 
structed from the coefficients b' m . 

Section l4~2l makes extensive use of equations <35M36t for 
each halo in our simulations to characterise statistically these two 
fields. 



3.3 From flux densities to the 5D source 

The description of the source term s e involves time dependent co- 
efficients c l ^, m ,{t). Their computation from the particles coordi- 
nates is quite straightforward and as shown in the previous sections, 
the different margins can be recovered through the manipulation of 
theses coefficients. Yet a projection of the source on an a priori ba- 
sis is a complex operation. Here this projection aims at describing 
a 5D space for which little is known. As shown in the following 
sections, the distribution of incidence angles is quite smooth, while 
the distribution of velocities appear to be easily parametrized by 
gaussians. The 5D basis presented in the current paper induces lit- 
tle bias, but it is very likely that a more compact basis exists and 
that the size of the expansions chosen can be reduced in the future. 

Because of this large amount of information contained in the 
source, it is not always convenient to relate coefficients or their 
correlations to physical quantities, like the mass flux or the flux 
density of e nergy. An alternate description of the source term was 
presented in lAubert et alj 120041) with the following ansatz: 



= (r,v,t) = ]TY m (ft) 



• (2t) 



-3/2 



dct|iz7 c 



(37) 



exp 



This representation of the source is by construction consistent with 
the first two velocity moments: 



while 



r, v = uj 



d 3 v( Vi 



fo- 



ci vvs e (r, v) = TOp V (r) 



)s e (r,v) = 



(38) 



zu pv (r) 2 
zu p (r) 

,(r). 



+ J2 Y m (fi)5(r - R 20 



»(*) 



(39) 



Obviously, the third moment is not fully recovered from the Ansatz 
given by equation 1381 . This example should be taken as an il- 
lustration and highlights the possibility of building a source term 
from its moments. It is not unique and more realistic expressions 
could be found, which satisfies higher moments of the source. Still, 
the successive measurements on the sphere of the flux density of 
mass zu p , momentum zu pv and velocity dispersion ■co paa allow 
a coherent description of the infall of matter. Unlike the coeffi- 
cients, these flux densities are easier to interpret since they describe 
physical quantities and are directly involved in specific dynamical 
processes (see table 13. 3t . Furthermore, these three flux densities 
are easily expressed in terms of coefficients c„™ m /(t), or more 
precisely in terms of a subset of the source's coefficients, imply- 
ing a smaller number of computations relative to a complete cal- 
culation of c a ™ m /(t). Finally these flux densities are particularly 
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suited to the regeneration of synthetic environments. As shown in 
appendix [F| synthetic spherical maps can be generated from the 
two-point correlations and cross-correlations of these fields. Such 
environments would be consistent with the measurements in sim- 
ulations and will allow us to easily embed simulated galaxies or 
halos in realistic environments as a function of time. 

The expression given in equation 1381 has one important 
drawback : it is not of the form of equation 1121 . i.e. it would 
require a reprojection over a linear expansion for a dynamical 
propagation. Nevertheless, its compacity makes it easier to com- 
pute than the full set of coefficients and the associated strategy 
would be 1) to measure the flux densities from the simulation, 2) 
build a source term from e.g. the Equation <38> and 3) project over 
an appropriate 5D basis when needed, i.e. when the source is used 
as an input to the analytic description of the haloes dynamics. 

The following sections will make intensive use of the coeffi- 
cients described by the equation 1351 and equation <36t . In particu- 
lar, it will shown how the manipulation of these coefficients allow 
to recover relevant physical quantities. In the current paper, only 
the first moment of the source, the flux density of mass vj p together 
with the external potential, will be fully assessed. The kinematical 
properties of the accreted material will in particular be investigated. 
The complete characterization of the c(t) coefficients is beyond 
the scope of the current paper and will be completed in paper III. 
The full measurements of these 1 1 fields required by Eq. l38l and the 
comparisons between the two expressions of the source will also be 
assessed in this future paper as well. Appendix lG 1 l and lG2l describe 
how the other moments, the flux density of momentum vj pv or the 
flux density of energy vj pa i may be recovered from the source ex- 
pansion. 
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Figure 5. An example of accretion history. Top : Azimuth-time diagram. 
Each point in the diagram represents one particle passing through the virial 
sphere at a given azimuth (y-axis) at a given time measured from the Big- 
Bang (x-axis). Bottom: the distribution of particles' crossing instants (in 
black). Time increases from left to right. The infalling (resp. outflowing) 
particles distribution is shown in red (resp. blue). 



Since 



&YY t , tm , = \/An5 t , 8 m , , and / dvv 2 g a (v) = ^ 2 +<7 2 .(41) 



It follows that 



(42) 



3.4 A template halo 

As an illustration, let us first apply the whole machinery to one typi- 
cal halo. At z=0, this 'template' halo has a mass Mof3.4-10 13 M Q , 
with a virial radius -R200 of 800 h^ 1 kpc. The corresponding cir- 
cular velocity is V c o = 600 km/s. Its accretion history is shown 
in figure [5] for z < 1. Each point on the azimuth- time diagram 
represents one particle of the simulation passing through the virial 
sphere at a given azimuth and at a given time. Temporal space has 
been sampled using 15 equally spaced bins between z=l and z=0 
(see bottom panel in figure [3}- For each time step, the expansion 
coefficients c™ m (t) are computed from equation 1351 . The Gaus- 
sian basis g a (v) involved 25 functions with means /i Q equally dis- 
tributed from v=0 to v=1.5 in V c o units and with an r.m.s a — 0.03. 
The harmonic expansions were carried up to £ — 50 in position 
space and up to £' = 25 in velocity space. 



3.4.1 Advected mass: angular space (vj p ) 

The template halo accretes an object at t s =l 1 Gyr (where t=14 Gyr 
stands for z=0) adding 7.5- 10 12 Mq to the system during a ~ 1 Gyr 
interval. The corresponding spherical flux density field, uj p (tt, t 3 ) 
is shown in Fig. {6}. It represents the distribution of accreted parti- 
cles as seen from a halo-centric point of view. The field uj p (fl, t 3 ) 
has been reconstructed from the coefficients (see equation i35\ ). It 
reads: 

w p {n,t) = J dTdvv 2 s e (v,n,T,t) = ^a m (t)F m (fi). (40) 



allowing us to recover zu p (fl, t s ). Also shown, the same field but 
computed this time using directly the a ngular distribution of parti- 
cles as described in lAubert et al. I l2004l) (see below). All the major 
features are well reproduced by the expansion coefficients, equa- 
tions I35i - I42l . Clearly, an object is 'falling' through the virial 
sphere. It is straightforward to obtain the angular power spectrum 
Cj 7p from the c™ m , coefficients via the definition of ai m in equa- 
tion J42J : 



c: 



i 



4tt 21 + 1 ^ 



(43) 



The angular power spectrum of uj p (fl,t s ), derived from the ex- 
pansion, equation <35> . is shown in Fig. 0. From the positions and 
velocities of particles it is also possible to evaluate vj p (ft,t 3 ) on 
an angular grid and recover the angular power spectrum 'directly' . 
The agreement between the two C^ p is good, though for the small- 
est scales {£ ^ 30), the power spectrum computed from the coeffi- 
cients is slightly larger than the one derived directly from the par- 
ticles. This may be explained by the fact that a grid sampling tend 
to smooth the actual zu p field. As a consequence the amplitude of 
small scales fluctuations is decreased, leading to a smaller C e p . A 
more complete discussion on harmonic convergence can be found 
in appendix [3] For a given £, the corresponding angular scale is 
n/£in radians. 

Note that the coefficients aoo are closely related to the accre- 
tion field averaged over all directions, $ M (t), defined by : 



(44) 
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flux density , ■w 



Flux, $ Motivation 



Mass pv r 
Angular momentum pv r r X v 
Kinetic energy pv r ai Oj 



dm I At heating & cooling, 
dL/dt warp, shape of halos, 
dE I dt virialized objects, 



Shear 
Vorticity 



pv r (dvj /dxi + dvi/dxj) dc/dt 



tidal field, 



pVr'V X V 



du/dt anisotropic accretion. 



Table 1. Description of the various flux densities. The first 10, together with the external potential are sufficient to characterize 
fully the environment as shown in Section l3~3l 
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Figure 7. The angular power spectrum, Ct/Ci, (equation 1431 ) of the dis- 
tribution of incoming matter, zu p (Si) (shown in figure |6j at t ~ 11 Gyr, 
for our template halo. For a given £, the corresponding angular scale \sir/£. 
The histogram corresponds to the power spectrum derived directly from the 
particles' angular distribution. The solid line is the power spectrum recon- 
structed from the coefficients. 



Figure 6. An example of a flux density reconstructed from the coefficients 
c™ m : the mass flux density, piy (fi). It represents the angular distribution 
of incoming mass as seen from a halo-centric point of view. Here t s ~ 1 1 
Gyr. Light regions correspond to strong infall while darker regions stand for 
low accretion and outflows. Top : the spherical field obtained directly from 
the spatial distribution of particles. Bottom: the reconstructed spherical field 
from the coefficients, equation 1351 . 



Projections over T2 and v gives the probability distribution of the 
incidence angle Ti, ^(Ti, t), defined as: 

#(ri,t) = / dr 2 dvv 2 pv r (r,v,t) , 



W5F£]tS, { i',o>0£ +<r 2 )Y e , i0 (T). (47) 



Measuring aoo (t) amounts to measuring the accretion flux density, 
i.e. the quantity of dark matter accreted per unit surface and per unit 
time. 



3.4.2 Advected mass: velocity space 

Integration over the sphere leads to the distribution of accreted mat- 
ter in velocity space : 



pv r (r,v,t) 



7r$>°, m , 5Q (i,)y m ,(r). 

a ,m' 



(45) 
(46) 



The impact parameter b of an incoming particle (measured in units 
of the virial radius) is related to Ti by 



-R200 



= sinfTi) , 



(48) 



therefore the probability distribution of impact parameters, #(&), is 
easily deduced from equation 1471 . At t ~ 11 Gyrs, the #(&) com- 
puted from the source coefficients is compared to that derived di- 
rectly from the velocities of particles in Fig. {8j. Note that for pure 
geometrical reasons small impact parameter b are less likely since 
there is only one trajectory passing through the center while there 
is a whole cone of trajectories with b 7^ 0. As a consequence er- 
rors are intrinsically larger for small values of b. The reconstruction 
from the source coefficients is clearly adequate. In this example, 
the high probability for infalling particles to have a small impact 
parameter (i> < 0.5) imply that velocities are strongly radial. The 
object 'dives' into the halo's potential well. 



Dynamical flows through Dark Matter Haloes II :One and Two points Statistics at the virial radius 1 1 



Projection over T\ leads to the probability distribution of par- 
ticles velocities, ip(v,t s ), as they pass through the virial sphere. 
The PDF tp(v, t) is defined as: 

(p(v, t) = v 2 J d£ldTs e {v, ft, T, t). (49) 

Here the v 2 weighting accounts for the fact that the probability dis- 
tribution of measuring a velocity, v within dv is of interest here. 
Using coefficients, it follows that: 

tp(v,t) =4ir^2v 2 g a (v)c° a , . (50) 

a 

The reconstructed velocity distribution is also shown in Fig. 
It reproduces well the actual velocity distribution. For this specific 
halo, the satellite is being accreted with a velocity of 0.75 V c o- 

The correlation between the incidence angle Ti and the ve- 
locity's amplitude v may be studied by integrating pv r (T,v,t) 
over T2 only. The distribution function, p(Ti,v), of particles in 
the (Ti, v) subspace is defined by: 

p(Ti,u) EE J dr 2 dvv 2 pv r (T,v,t) , 

= 2jrV4jF5^c£, {i ,,o } s«(«)y//,o(ri,0). (51) 

a,i' 

Given the relation equation 1481 . the correlation p(b, v) between 
the impact parameter and the velocity's amplitude is easily ob- 
tained. The p(b, v) distribution is shown in Fig. J8}. Again, note 
that p(b,v) represents an excess probability of finding an impact 
parameter b (with a velocity v) compared to isotropy. In this specific 
example, no real correlation may be found between the two quanti- 
ties. Finally, the integration of p(b, v), tp(v, t) and $(Ti) over their 
respective space leads to the same quantity, namely the integrated 
flux<I> M (t). 

3.4.3 External potential 

The final field needed on the virial sphere is the external tidal field 
created by the dark matter distribution around the halo. 

Using equation 1361 . the external potential tp e (ft, t) is easily 
computed from the positions of external particles, having restricted 
the sampling to particles within a 4 Mpc (physical) sphere centered 
on the halo. The position of external particles are linearly inter- 
polated at a given measurement time. The b m coefficients for the 
template halo are computed at t„ ~ 7 Gyrs (measured from the 
Big Bang). The reconstructed ip e (Cl,t) field is shown in Fig. |9j 
along with the modulus of the advected mass |piv(f2)|. The two 
reconstructions were restricted to harmonics i ^ 20. 

The two spherical fields show the same main features. How- 
ever, almost no small scales features is seen in the map of the ex- 
ternal potential even though they have the same resolution. Since 
the gravitational potential is known to be smoother than the asso- 
ciated density and is dominated by the global tidal field, it is not 
surprising that ip e (fl, t) appears smoother than the advected mass 
field |c7p(fl)|. 

The potential's angular power spectrum may also be computed 
by replacing ae. m by be m in equation 1431 (see Fig. J 101 ). The power 
spectrum of the potential, C e p , sharply decreases with I, while 
Cj 7p has a gentler slope. Large scales are clearly more important 
for the potential than for the advected mass. Furthermore, Cf sys- 
tematically peaks for even i values, reflecting the 'even' symmetry 
of the potential measured on the sphere. 
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Figure 8. Top: Excess probability distribution of impact parameter b, t9(t>), 
derived from the c^ m , (t) source coefficients (equation 1351 ) of our tem- 
plate halo at t s ~ 11 Gyr (line). The histogram corresponds to the same 
distribution derived directly from the particles' positions and velocities. Er- 
ror bars stand for 3<r errors. The impact parameters are given in units of 
i?200- The Y-axis unit is 5 ■ 10 9 MQ/kpc 2 /Myr. Infall is mainly radial. 
Middle: The velocity distribution of particles, ip(v), accreted at t s ~ 11 
Gyr, for our template halo. Velocities are expressed in terms of the circular 
velocity at z=0. The Y-axis unit is 5 ■ 10 9 MQ/kpc 2 /Myr. The histogram 
corresponds to the velocity distribution obtained directly from the particles 
velocities. The solid line is the reconstructed distribution from the source 
coefficients. Bottom: the probability distribution, p(b, v). of particles in the 
b — v subspace. Units are the same as above. Red/blue stand for high/low 
densities. No correlation is found between b and v for this specific example. 
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4 SIMULATION SAMPLE & STATISTICAL BIASES 

Section [3T4l details the measurement strategy for a given typical 
halo. It is now possible to reproduce the above measurements for all 
the halos of the simulation sample. Let us first describe in turn the 
construction of our sample, and the corresponding biases, which 
constrain our ability to convert a large set of simulations into the 
statistics of the source. 



4.1 Simulations 




Figure 9. A comparison between the external potential, »/) e (f2), and the 
modulus of the flux density of matter, \pv r (ii)\. The measurement is made 
at t ~ 7 Gyrs (measured from the Big Bang) on our template halo. The two 
fields were respectively reconstructed from c™ , and fe m coefficients with 
^max Ss 20. Even though the two fields are similar and exhibits a strong 
quadrupolar component, i/! e (f2) is smoother than pv r (S~l). It is expected 
that the corresponding expansion coefficients be statistically correlated. 




harmonic L 

Figure 10. A comparison between the angular power spectrum of pv r (fl) 
and i/) e (f2) for our template halo (the two fields are shown in Fig. The 
two power spectra C'i are normalized by C%, i.e the quadrupole contribu- 
tion. The slope of the potential's power spectrum is clearly stronger. Large 
scales (i.e. small I values) dominate the angular distribution of i/> e (f2), as 
expected. 



In order to achieve a sufficient sample and ensure a convergence 
of the measure ments, a set of ~ 5 00 simulations was produced 
as discussed in I Aubert et alj l2004t) . Each of them consists of a 
50 h^ 1 Mpc 3 box containing 128 3 particles. The mass resolution 
is 5 • 10 9 M Q . A ACDM cosmogony (fi m = 0.3, fi A = 0.7, 
h = 0.7 and cr 8 = 0.928) is implemented with different initial 
conditions. These in itial conditions were produced with G RAFIC 
iBertschingeilfeOOll) ) where a BBKS iBardeen et at] i 19861) 1 trans- 
fer function was chosen to compute the initial power spectrum. 
The initial conditions were used as inputs to the parallel version 
of the tree code GADGET ISpringel et ail feOOll) ). The softening 
length was set to 19 ft -1 kp c. 3 The halo detect i on wa s performed 
using the halo finder HOP jEisenstein & Hutl jl998l) 1. The den- 
sity thresholds suggested by the authors (<5 ou ter — 80, <5 sa ddic = 
2.5<5 ou tcr, <5p Ca k = 3.<5 ou tcr) were used. 



4.2 Selection criteria 

As shown in lAubert et alj <2004 the completion range in mass of 
the simulations spans from 3 • 1O 12 M to 3 • 10 14 M o . Since the 
emphasis is on L* galaxies, the survey is focused mainly on galac- 
tic halos and light clusters, only haloes with a mass smaller than 
1O 14 M at z=0 were considered. The interest is for halos already 
'formed', i.e. which will not experience major fusions anymore. To 
satisfy these requirements the focus is on the last 8 Gyrs (redshifts 
z < 1 in a ACDM cosmogony). Since the history of a given halo is 
followed by finding its most massive progenitor, it is required not 
to accrete more than half its mass in a two-body fusion. As a final 
safeguard a halo is rejected if it accretes more than 5 • 10 12 Mq be- 
tween two timesteps (i.e. per 500 Myrs, see next subsection). This 
mass corresponds approximatively to the smallest haloes consid- 
ered at z = 0. The final range of mass of haloes which satisfy these 
criteria is ~ 5 • 10 12 Mq — 10 14 Mq, the fraction of rejected haloes 
being ~ 20%. Clearly, such a priori selection criteria will modify 
the distributions of measured values and the related biases may be 
difficult to predict. For instance, the Fig. lilt shows the scatter plot 
of the contribution of — 4) = 45 degrees fluctuations to w p 
field vs aoo, i-e. the accretion rate. It appears from this plot that 
modifying the treshold for the accretion will modify the average 
angular scale of w p in a non trivial way. Since only a small frac- 
tion of haloes is rejected, the biases are expected to be moderate, 
but as for now their impact cannot be estimated accurately on the 
average source or its moments. 

Let us emphasize that the above selection criteria should be 
added to those corresponding to the simulations themselves. Aside 
from the fact that a 50Mpc 3 /i _1 box size implies a limited range 



3 A second set of simulations with a resolution increase by a factor 2 3 
(resp. 2 6 ) was carried in order to investigate the convergence of some mea- 



surements (see Sec. 16. 2. 21. 
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Figure 11. Scatter plot of C™ p versus a^ Q measured at lookback times t= 
7.8 (red), 5.6 (green), 4. (yellow) and 2.9 (blue) Gyrs. The quantity aoo 
scales as the average accretion rate of the haloes while C 4 p scales as the 
contribution of I = 4 structures in the flux density of mass measured on the 
sphere. This plot illustrates how a threshold on the accretion rate affects in 
a non trivial way the typical clustering measured for zu p . In particular, one 
should note how C 4 p remains constant at recent times for low accretion 
rates. 



of mass, the universe described in these simulations is more ho- 
mogeneous than it should be, since each box must satisfy a given 
mean density. In other words, the probability of rare events is re- 
duced. This effect should not influence the number of haloes with 
high accretion rate, since strong accretions are rejected a priori. On 
the other hand, it should influence the number of objects which ex- 
perience low accretion history, which are probably less numerous 
in our simulations than in larger simulated volumes since voids are 
less likely. Furthermore, the intrinsic mass resolution sets a mini- 
mum accretion rate equals to one particle mass (5 • 1O 9 M0) per 
time interval. One could imagine an object with a mass smaller 
than the particle's mass which would not be included in the simu- 
lation at the current resolution. Furthermore, an object with a mass 
equals to a few times the minimum mass is be considered as diffuse 
accretion. Finally this mass resolution is related to the spatial res- 
olution, which limits intrinsically the angular description of fluxes 
on the virial sphere. For a given type of simulation, all these effects 
cannot be avoided and reduce the representability of the following 
measurements. 

In short, this strategy involves a bias in mass, in redshift, in 
resolution and in strength of merging event. However these biases 
should only influence somewhat extreme realisations (related to 
e.g. very low accretion or equal mass mergers) of the source or 
the external potential and since the focus is on the typical scales, 
presumably related to moderate interactions, they should hopefully 
not significantly affect the measurements. 



4.3 Reduction procedure 

In the following discussion most of the distances (resp. velocities) 
will be expressed as functions of the virial radius -R200 (resp. the 
circular velocity V c ) measured at z = 0. These quantities are re- 
lated to the halo's virial mass by 
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V c =,f^. (52) 
V W200 

Here A/200 = M(r < i?20o)- The mass-dependence of -R200 an d 
V c are given in Fig. 1 1 2i and may be fitted by: 

#200 = 537Af 1/3 , V c = 400M 1/3 , (53) 

where -R200 is expressed in (/i _1 kpc), V c in km/s and M in units 
of 1O 13 M0. Here M stands for the total mass of the halo, returned 
by the halo finder HOP. In Equation <53> . -R200 and Vc appear to be 
strongly correlated to the final masses of haloes, the few outliers be- 
ing related to external subhalos or peculiar halo geometries. Since 
the selection criteria are quite restrictive, most of the halos experi- 
ence the same relatively quiet history of accretion and account for 
the lack of scatter. 

The simulations in that redshift interval involved 15 snapshots 
sampled with A(logz) = est for z ^ 1 down to z = 0.1 plus a 
snapshot at z — 0. The gap between the last snapshot and the sec- 
ond to the last is nearly 1.4 Gyr. As a consequence, the assumption 
of ballistic trajectories is not valid anymore (see appendix) . Simu- 
lations were re-sampled in 15 bins distributed regularly in time (i.e 
not in redshift) using the procedure described in Section l3~2l the 
corresponding time step is ~ 500 Myrs. To take in account the last 
gap, results obtained from the last three 'new' bins (which cover 
the last 1.4 Gyr) were averaged into a single bin centered on 0.8 
Gyr. 

The source coefficients, c™ m /(i), were computed following 
the procedure described in section l3~4l Maximum harmonic orders 
were set to i? max = 50 for the position-angular description. For a 
typical halo with ifeoo = 500 kpc, l max = 50 corresponds to a 
spatial scale of 30 kpc, i.e. equals to 1.5 times the spatial resolution 
of the simulation. The harmonic description of the velocities angu- 
lar dependence is restricted to l' max = 15. The velocity amplitude 
is projected on a gaussian basis which involves 25 functions regu- 
larly spaced from v/v c = to v/v c = 1.5 with an r.m.s of 0.03. 
These parameters allow a satisfying reproduction of distributions 
computed from particles. 

The external potential coefficients, b m (t), were computed fol- 
lowing the procedure described in section l3~4l Only particles within 
a 4 Mpc physical sphere centered on the halo are taken in account. 
Maximum harmonic orders were set to t max = 20. 

A set of 100 simulations have been fully reduced allowing us 
to compute c™ m , (t) and b m (t) for 15000 halos. Since a well de- 
fined, (if only biased) sample of histories of haloes was constructed 
in our simulations, it may be projected on our basis, to compute 
the external potential and the flux density of mass, following Sec- 
tion l3.4l Let us now characterize the corresponding coefficients, via 
one point (Section|5J and two points (Section|6j statistics. 



5 ONE POINT STATISTICS 

In this section, let us first describe the evolution and the statisti- 
cal distributions of the global properties (i.e. integrated over the 
sphere) of the source and the potential. Let us discuss the evolution 
of the mean potential, of the mass flux and the kinematical 

properties of s e via the velocity distribution (f>(v) and the impact 
parameter distribution $(&). 

Let us then describe in turn the statistical PDF, mean and vari- 
ance of the integrated fluxes, their corresponding flux densities, and 
finally their mean kinematical features, following the steps of sec- 
tionEH 
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Figure 12. Virial radii (-R200) and circular velocities (V c ) as functions of 
halos final masses. The quantities have been measured at redshift z = 0. 
Scaling relations between R200 or V c and the final mass are also given. 



5.1 Mean External potential 

The mean external potential on the sphere is actually somewhat 
meaningless but is being used as a normalization value for potential 
fluctuations (see section loTV Because of isotropy, the mean poten- 
tial is seen as a monopole and only the 600 (t) coefficient is statisti- 
cally different from zero. Furthermore, following equation J24> . the 
three dimensional potential component induced by the monopole 
is a constant potential throughout the sphere volume. As a conse- 
quence, it influences the halo's dynamics only through its temporal 
variation. 

The time evolution of the (600)} coefficient is given in 
Fig. 1 1 3i . The &oo(i) distribution exhibits a tail due to large —600 
values and is better fitted with a log-normal distribution than with 
a Normal distribution (see Fig. IC1L Hence ((boo)} stands for the 
most probable value of the log-normal fitting distribution. 4 The 
evolution shown in Fig. 1131 reflects the measurement procedure. 
Given that the potential is computed from all the particles contained 
within a fixed physical volume, the overall expansion implies that 
particles tend to exit the measurement volume with time. In other 
words, the average density in the measurement volume decreases 
with time. This effect leads naturally to the decline of the average 
potential within the virial sphere due to external material. 



5.2 Mass flux : $ M (i) 

At each time-step, the ao.o distribution is fitted by a Gaussian func- 
tion with mean (eto,o) (see also Fig. <C3» . This gaussian hypothesis 
is clearly verified for low redshifts while strong accretions events 
give rise to a tail in the ao,o distribution at high z. At these epochs 
(lookback time t > 7 Gyr) the gaussian fit tends to slightly over- 
estimate the mode position. Yet the gaussian hypothesis remains a 
good approximation of the distribution while the time evolution of 
the gaussian mean value (ao,o)(i) represents well the evolution of 
the mode of ao o- 



4 Since the measured distribution is quite peaked, fits made with a Normal 
distribution (not shown here) return a very similar time evolution of the 
mode position. 
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Figure 13. The time evolution of the monopole component of the external 
potential, &oo(i). The time evolution is fitted by a second order polynomial 
-b 00 (t) = 35948 *t 2 - 61480.7 * t + 793067. Look back time t is 
expressed in Gyrs while &00 coefficients are expressed in units of GM/R. 
M is expressed in 10 10 Mq, R in kpc/i -1 and G=43007 in internal units. 



The time evolution of the average flux of matter through the 
sphere, <t M (t) = is directly derived from the evolution of the 
monopole (see equation <44> 1 and is shown in Fig. J 1 4-1 . It can be 
fitted by: 



$ M (t) = <7cv)(i) = -0.8H 3 + lQ.lt 2 



19.3* + 17.57, (54) 



where <& M (i) is in units of M0/Myr/kpc 2 and look-back time t 
is expressed in Gyr. As expected, the average quantity of material 
accreted by halos decreases with time. For z < 1, a large frac- 
tion of the objects of interest in are already 'formed' and only gain 
matter through the accretion of small objects or diffuse material. 
In a hierarchical scenario, such source of matter become scarcer, 
inducing a decrease in the accretion rate. Furthermore, recall that 
$ M (t) is measured as a net flux, i.e. the outflowing material may 
cancel a fraction of the infalling flux. Therefore, the decrease with 
time may also be the consequence of an increasing contribution of 
outflows: the measurement radius i?2oo(z = 0) becomes the actual 
virial radius of the halo as time goes by, i.e. the radius where the 
inner material is re-processed and where outflows are susceptible 
to be detected. 

As a check the average mass accretion history (MAH here- 
after) of the haloes was computed. The MAH ^(t) is defined as: 



*(*) 



M(t) 



(55) 



M(z = 0) 

where M(z) is the halo's mass at a given instan t . Using the ex- 
tended Press-Schechter formalism. lvan den Boschl l2002ah showed 
that halos have an universal MAH, fitted by the following formula: 



log(*(A/(z = 0),t)) = -0.301 



Ml + z) 
log(l + z/) 



(56) 



where z / and v are two parameters which depend on the considered 
class of mass on ly. These two paramete rs are found to be correlated 
and for instance [Wechsler et alji2002T) found a similar relation us- 
ing a single parameter. For each halo, its mass evolution M (t) was 
computed from its final mass M{z = 0) and its integrated flux of 
matter $ M (t): 
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Figure 14. Top: Time evolution of the average flux density of matter 
through the virial sphere, ( < I >M (t)) = (za p )(t) (symbols). Bars stand for 
3<t eiTors. Here <J> M (t) is computed directly from ooo coefficients follow- 
ing equation 1571 . Its time evolution is fitted by a 3rd order polynomiaKso/irf 
line). Bottom: the mass accretion history log M(z)/M(z = 0) for three 
different class of masses. Masses are expressed in solar masses. Symbol 
represent the median value of log M(z)/M(z = 0) within each classes . 
Lines represent the fitting function suggested by Ivan den Bosch] <2002j) . 
Even though the global behavior is reproduced by the fitting functions, the 
measured accretion rate is systematically sm aller. This discrepancy has al- 
ready been noted bv lvan den BoschH2002al) . 



M(t) = M(t = 0) - A-kRI 



dt$ M (t). 



(57) 



From equation J55I ^(t) was computed for each halo. The median 
value of fyjt ) for three classes of m ass was compared to the fit 
suggested bv lvan den Boschl l2002at) (see Fig. <14» . For the three 
classes, the agreement with the fitting formula is qualitatively satis- 
fying: the three measurements evolve in the expected manner while 
their relative positions are the same as the relative positions of the 
three fits. However our measurements are quantitatively inconsis- 
tent with the three curves. At low redshift, ^(t) is systematically 
larger than the expected value (i.e. the accretion rate is smaller). 
The median mass at z = 1 is well recovered even though the two 
method disagree slightly quantitatively. In other words, our mea- 
surements overestimates the accretion at high redshift and underes- 



timates accretion at low redshift. This may be related to the mea- 
surement procedure through the sphere : at higher redshift, accreted 
material is assumed to be added to the biggest progenitor, even 
though it has not yet reached the central object and its mass is over- 
estimated. Still, this material end up in the most massive progenitor 
and the final mass is recovered. Note that since specific selection 
criteria were applied, these halos may not be completely represen- 
tative of the whole halos population. Finally recall that the median 
value of ^(t) was represented here because of strong outliers while 
the fitting formula is given for the average MAH (extracted from 
merger trees). A similar discrepancy between the extended Press- 
Schechter theory and the results obt ained from numerical s imu- 
lations had already b een noticed byjrar^ei^oschjj]2002a|) and 
IWechsler etaD <2002l) . In particular, van den Bosch 1 2002a) found 
that the Press-Schechter models tend to underestimates the forma- 
tion time haloes compared to simulations. Clearly, our measure- 
ments seem to confirm this discrepancy. Since the global behavior 
of MAHs is recovered and since the median mass at z — 1 is recov- 
ered, it is concluded that the measure of $ M (t) through the virial 
sphere reproduces the accretion history of halos. 



5.3 Mean kinematics 

Let us now turn to the kinematical properties of the flow, while 
averaging the source over the virial sphere. 



5.3.1 probability distribution of the modulus of velocities 

Given the source coefficients c™ m , , the average velocity distribu- 
tion (<p(v,t)} (defined by equation <49» is easily computed since 
it only involves (c° j0 (i)). The ensemble average of (Ca,o) and the 
related ensemble dispersion cr(Cao)(t) = (( c a,o ~ ( c a,o)) 2 ) aK 
derived by fitting the c° a (t) distribution by a Gaussian function. 
From these two quantities, it follows: 



(ip(v,t)) = 4irv 2 ^ gq(«)(c° [0 ) , 



and 



ff [ v »(«,t)] = W /X>MM<o) 2 , 



(58) 



(59) 



which are respectively the ensemble average and r.m.s. of the veloc- 
ity distribution. The time evolution of {<p(v, t)) is given in 

Fig.O 

and 1161 . Errors on {<p(v, t)) are computed as 



A [<*>(«,*))] = 3 



OS 



(60) 



At 'early times' (t > 5 Gyr), the distribution is unimodal with 
a maximum around 0.7V c (z — 0). No outflows can be detected at 
any velocity and the infalling dark matter dominates. At later times, 
(<p(v,t)) drops below zero for velocities around 0.4V c (z = 0). 
Outflows dominate at 'low' velocities. Meanwhile the amplitude 
of the previous peak decreases and shift to higher velocities. The 
fraction of infall relative to the total amount of material passing 
through the sphere drops from 1. to 0.6 between t = 8 Gyr and 
t = 0.8 Gyr. 

This behavior is likely to be due to our measurement at a. fixed 
radius, R.2oo(z = 0). At 'early times', this measurement radius is 
bigger than the actual virial radius of halos. Thus no sign of 'viri- 
alization' (outflows consecutive to accretion) is detected. Later, the 
actual ifeoo gets closer to the measurement radius. Outflows pass 
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Figure 15. The time evolution of the average velocity distribution, 
(ip(y,t)), defined by equation 1491 . for z < 1 (symbols). Ages are ex- 
pressed as look-back times (i.e. t = for z = 0). Velocities are 
given relative to the halo's circular velocity at z = 0. Y-axis unit is 
5 ■ 10 9 MQ/kpc 2 /Myr/V c . Error bars stand for 3 — a errors. Here <fi(v) 
is fitted by the sum of two Gaussian with opposite signs (solid line). Each 
Gaussian's contribution is also shown (dashed lines). 



through the measurement radius as a sign of internal dynamical 
reorganization. The fact that accretion intrinsically decreases with 
time would provide another explanation of this trend. This decrease 
can actually be traced in Fig. <15> . 1 1 61 and 10 1 1 . 

The global behavior of {<p(v, t)) can be modeled by summing 
two Gaussians representing the infalling and outflowing compo- 
nents : 



MM)) 



g.,3(t) 



<?i,2(i)v27T 

9o,s(t) 



exp 



+ 



<7o,2(i)v27r 



exp 



(61) 



Subscripts i and o stand for infall and outflow. Note that (t) ^ 
and q ,s(t) ^ 0. The coefficients time evolution is given in 
Fig. JDlk where t should be expressed in Gyr and <p(v, t) unit is 
5 • 10 9 M Q /kpc 2 /Myr/ V c . Examples of fits are shown as solid lines 
in Fig. <15t . Note that all the six coefficients evolve roughly linearly 
with time (see Fig. iDli ). Their linear fitting parameters are given 
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Figure 16. The time evolution of the average velocity distribution in the t- 
V/V c plane. Red symbols stand for positive values of the distribution (i.e. 
infall) while blue symbols stand for negative ones (i.e. outflows). Each of 
these components is fitted by a Gaussian function in the V/V c space. The 
time evolutions of the mean of the Gaussians are given by the two lines 
(solid for infall, dashed for outflows). The Gaussians r.m.s. are also shown 
as bars. 



in table lDll Using equation 16 1 1 and the linear parameterization of 
the Gaussians coefficients, {f(v, t)) is reproduced accurately. The 
only restriction concerns the negative Gaussian's amplitude (q ,3) 
which should not be greater than 0. Since this condition is not nat- 
urally satisfied by a linear fit, it should be set by hand. 

The evolution of the relative positions of two Gaussians is 
given in Fig. 1161 . For t > 5Gyr, it is consistent with the 'no 
outflow' hypothesis, the amplitude of the negative Gaussian being 
close to zero at this epoch. Both Gaussians mean values seem to 
drift to higher velocities as a function of time. Even though the rel- 
ative velocity of accreted material is determined by the initial con- 
ditions (namely large scale clustering), the velocity of an infalling 
satellite should partly reflect the properties of the accreting body. A 
dense massive halo will not accrete like a fluffy light one. In other 
words, the velocity of infalling material should reflect the actual 
circular velocity of the accreting body. As a consequence it is ex- 
pected that accretion velocity drifts with time towards V c (z — 0). 

Furthermore, the mean values of both Gaussians evolves 
roughly linearly over the whole time range with comparable rate 
of change (see Fig. i Q 1 1 and the following discussion). As a con- 
sequence, their relative positions remains roughly constant (see 
Fig. <16» . This indicates that these two components may be phys- 
ically related, outflo ws being the consequence of a past accretion. 
iMamon et alj 120041) mention the existence of a backsplash popu- 
lation, rebounding through the virial radius and t his po pulation is 
known to have a different velocity (e.g. lGill et alll2004l) ). The out- 
flows detected via our description of the source is consistent with 
this backsplash component. The difference in velocity may be ex- 
plained if outflows are representative of an earlier accretion with a 
velocity typical of earlier times. Also past accreted material is in- 
fluenced by the halo's internal dynamics. Its velocity distribution 
would be 'reprocessed' (e.g. via dynamical friction, tidal stripping 
or phase mixing) to lower velocities as the material exits through 
the measurement radius. 

However, recall that the distribution shown here is a 'net' dis- 
tribution. In other words, it is quite plausible that an outflowing 
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Figure 17. Top: the net velocities distribution (histogram) measured from 
300 haloes at t=3.4 Gyr. This distribution is representative of the distribu- 
tion computed from coefficients and averaged over 15 000 haloes. Velocities 
are given in circular velocity unit while y-axis units are arbitrary. The two 
components are fitted by two gaussians (dashed lines). Bottom: the sepa- 
rate velocities distribution of accretion (red histogram) and outflows (black 
histogram). The dashed curves represent the difference between these two 
distributions and their respective fits shown above. It results in two residual 
distributions, centered on the same velocity and displaying nearly the same 
shape. These two residual distributions describe the material which already 
experienced one passage through the virial sphere. 



component may be completely canceled by an infalling component 
which has an exact opposite distribution. This effect is illustrated 
by Fig. < 171 where the velocity distribution of infall and outflows 
are being shown separately. This distribution has been computed 
from 300 haloes at t=2.3 Gyr. This distribution is quite representa- 
tive of the average ones, except a few high velocities events which 
skew the distribution of the infalling component and which are in- 
duced by outliers. If the gaussians fits are removed from these two 
separate distributions, two almost identical distributions appear for 
the two components, centered on V/V c ~ 0.6. These two identical 
distributions are related to the virialised component of infall, which 
already interacted with the inner region of the halo. The overall 
shape of these two distributions may provide insights on the typical 
dynamical state in the haloes' inner regions. 



5.3.2 Impact parameters and incidence angles 

The average distribution of incidence angle, ($(ri,t)} has been 
computed following the same procedure described above for 
(tp(v,t)). Defining c e /(t) as: 

Cf/(<) = 27p v / 47r^c ai0> {y >0 }(/4 + a 2 ), (62) 

a 

yields 

<0(ri,t)) = 5^Y</,o(r)(a e /(t)> > (63) 

v 

and 

<r(0(ri,t)) = ^^Y^ (T)Mc^))2, (64) 

where (i?(ri,t)} and <j($(Ti,t) are derived by fitting the c^{t) 
distribution by a Gaussian function. Errors on (#(ri, t)) are com- 
puted similarly to errors on {<p(v, t)) (see equation 1601 ) . The time 
evolution of (i?(ri,i)) is shown in Fig. J 1 81 . Since the impact pa- 
rameter and the incidence angle are simply related by b/ 'R200 = 
sin(Fi), (#(&, t)) is also easily computed. 

The infall (Ti > 7r/2 or the upper branch in (#(&, t)) dia- 
grams) is clearly mostly radial. The infalling part of the distribu- 
tion peaks for Ti ~ n instead of having a uniform behavior and 
this trend can be observed for all redshifts below 1. The distribu- 
tion slightly widens with but remains skewed toward large values of 
Ti. In the (#(&, t)) representation, the higher branch becomes flat- 
ter with time. The outflows ( Ti < 7t/2 or lower branch in ($(&, t)) 
diagrams) are mainly undetectable at early times, as mentioned ear- 
lier. As time increases, the outflow contribution becomes stronger 
and radial orbits (Ti ~ 0) also appear to be dominant. However 
the behavior of the 'outflowing' part of the {$(Ti, t)) distribution 
is almost linear and does not peak. Tangential orbits cannot be ne- 
glected for this component. 

The evolution of (#(Fi,i)) can be fitted by the following 
parametrisation: 

<tf(ri,t)> = J° - exp (- ( i) ~ rf ) +P2(t)r 1+P3 (t), (65) 

where po = 2 • 10~ 6 in our units. The 'infalling' part is modeled 
as a Gaussian while the 'outflowing' part is fitted linearly. The time 
evolution of the three parameters Pk(t), k = 1, 2, 3 can be fitted 
by a linear evolution and the related linear parameters are given in 
table ID2l The evolution of pi(t) confirms that the 'infalling' part 
of the distribution, {d(Ti, t)) widens with time. 

This result implies that material experiences a circularisation 
as it interacts with the halo. Consequently, orbits are more tangen- 
tial as particles exits and re-enter th e halo's sphere. S uch an effect 
has already been measured by e.g. iGill et all 120041) . Dynamical 
friction would provide a natural explanation f or this evolution of 
th e orbits, but this argument is refuted by e .g. IColpi et ail 1 19991) 
or lHashimoto et all 1200 3l) . IGill et all 120041) mention the secular 
evolution of halos to explain this circularisation: the time evolution 
of the potential well induced by the halo would affect the orbits of 
infalling material and satellites. Other processes such as tidal strip- 
ping or satellite-satellite interactions may also modify the orbital 
parameters of dark matter fluxes. Clearly, the interactions between 
the infall and the halo drive this circularisation but the detailed pro- 
cess still has to be understood. 

This dynamical circularisation could also explain why the 
'outflowing' part of the (i?(ri)) (or ($(&))) is flatter than the in- 
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falling one: by definition this component interacted with the halo in 
the past, unlike most of the infall. Finally, the Ti or b representation 
explicitly separates infall and outflows. It implies that virialized 
particles which pass through the sphere do not 'cancel' each other 
and do contributes to the distributions. Such a 'relaxed' material is 
likely to have a non zero tangential motion, flattening the distribu- 
tions as its contribution becomes important. Since the actual size 
of the halo gets closer to the measurement radius as time advances, 
this component contributes more with time and its flattening effect 
on the incidence angle (or impact parameter) distributions should 
increase as well. 

Fig. 1201 presents the correlation between the velocity ampli- 
tude v and the impact parameter b, at four different instants and for 
both the infall and the outflow component. Considering these two 
components separately, no correlation can be found: the incidence 
does not depend on the amplitude on first approximation. The only 
noticeable result comes from the fact that accreted material has sys- 
tematically a higher velocity than outflows, which confirms the re- 
sults obtained from the distribution of velocities only. Again, this 
effect is related to the separate origin of these two fluxes, accretion, 
being dominated by newly accreted material, and outflows, which 
were processed by the inner dynamics of haloes. 



6 TWO POINT STATISTICS 

Let us now focus on second order statistics, through the correlations 
on the virial sphere. The two point correlations are assessed through 
the angular power spectrum and the angulo-temporal correlation 
function for both the external potential, tj) e and the first moment of 
the source term, i.e. the flux density of mass, vj p . 

6.1 External potential 

6.1.1 Angular power spectrum 

The potential's angular power s pectrum Of is computed for each 
halo from the be m coefficients <Aubertetal]j2004l) ): 

b t , m = V4tt(^\ - 5«o-^t), (66) 
(ooo) (ooo) 

related to the potential contrast : 

6 m (SI) = - ^ = y b t , n Y t , m (n). (67) 

^ e > fx 

The probability distribution of Of (t) was weighted as described 
in appendix^ For each time step and for each harmonic I, the Of 
was fitted by a log-normal distribution (see Fig. lC2t . Let us define 
{Cf }} (t) as the mode of the fitting distribution. The time evolu- 
tion of the external potential's power spectrum is shown in Fig. 1211 
Globally the power spectrum is dominated by large scales and is 
quite insensitive to time evolution. However, two regime may be 
distinguished. For low order harmonics, {Cf } (t) remains mostly 
constant. For smaller scales (I > 5), {Cf }(t) increases along 
time. As a consequence the power spectrum's amplitude does not 
change but its shape evolves while smaller scales become more im- 
portant relative to larger scales. 

These two regimes reflect the twofold nature of tidal interac- 
tions of a halo with its environment. Small angular variations of 
the potential relate to small spatial scales and presumably track the 
presence of objects which are getting closer or going through the 
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Figure 18. The time evolution (symbols) of the distribution of the av- 
erage incidence angle Ti (^(Ti)), defined by equation 1471 . Ages are 
expressed as look-back time. Bars stand for 3a errors. Y-axis unit is 
5 ■ 10 9 M Q /kpc 2 /Myr. Outflows are counted negatively, leading to nega- 
tive values of (i9(ri)) for Ti < ir/2. The result of the model described in 
equation 1651 is also shown (red line). 

virial sphere. Since small scales contribution increases, it suggests 
that these objects tend to get smaller with time. It would be con- 
sistent with the global decrease of the accretion rate, as long as the 
merger rate does note increase strongly during this epoch. How- 
ever, the rise of small scales may also be related to an increasing 
contribution of weak and poorly resolved accretion events. In such 
a case, the isolated particles' contribution to the potential should be 
measured. This possibility is investigated in section l6~2l 

Meanwhile, large scale potential's fluctuations (£ ^ 4) may 
reflect the 'cosmic tidal field' resulting from the distribution of 
matter around the halo on scales larger than the halo's radius. The 
amplitude of such a tidal field should remain fairly constant, as in- 
deed measured. Furthermore, the peripheral distribution of matter 
is not spherically distributed but is rather elongated along some di- 
rection: haloes tend to be triaxial with their ellipsoid aligned with 
the surrounding distribution of satellites. The intersect of an elon- 
gated distribution of matter with the virial sphere would induce a 
quadrupolar component, as detected in our measurements. These 
two effects cannot be easily disentangled, since they actually are 
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Figure 19. The time evolution (symbols) of the impact parameter b distri- 
bution t)), defined by equation 1471 . Ages are expressed as look-back 
time. Bars stand for 3cr errors. Y-axis unit is 5 • 10 9 MQ/kpc 2 /Myr. The 
lower (resp. higher) branch is the t)) distribution for outflows (resp. 
infall). The result of the model described in equation 1651 is also shown (red 
line). 



two sides of the same effect. Large scale distribution of matter is 
responsible of both the 'cosmic tidal field' and the halo triaxial- 
ity (via the distribution of satellites). In other words, it is not clear 
whether the large scale behavior of (Cf } (t) reflects the tidal field 
or the halo's reaction to this tidal field. 
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Figure 20. Top: the distribution, (p(b, v)), of particles in the (6, v) sub- 
space at lookback time t = 1.8, 3.4, 5.1 and 6.7 Gyr; the infall (contour 
plot) and outflow (density plot) are represented separately. Beyond the bi- 
modal feature, no residual correlation appears. 
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Figure 21. The angular power spectrum of the external gravitational po- 
tential ?/) e (f2,t). Symbols represent the mode of the Cf distribution for 
each harmonic £ and each time step. Times are lookback times. Bars stand 
for 3 — a errors on the mode value. The large scales contribution remains 
constant with time while the small scales contribution smoothly increases 
with time. The bump for the (C2) component indicates a strong quadrupo- 
lar configuration for ip e (ft, t). 



6.1.2 Angulo-temporal correlation 

6.1.2.1 Correlations and coherence time To investigate further 
these two regime of tidal interactions, let us compute the angulo- 
temporal correlation function of the external potential contrast, de- 
fined as 

{w e (9, t, t + At)} = {8 m (O, t)S m (SI + Atl,t + At)}, (68) 
which is related to Tf (t, t + At) coefficients by: 

w e (6,t,t + At) = ^Tf (t,t + At)(2£ + l)P e (cos(6)), (69) 

e 

where 



Tf{t,t + At) = l-^—^2bi m (t)b* tm (t + At). (70) 

m 

Here, 9 stand for the angular distance between two points on the 
sphere located at and ft + ACl. The Tf (t, t + At) coefficients 
were computed for each halo and each pair of time step, for each 
harmonic. The Tf (t, t + At) distributions were fitted by a log 
normal distribution and {Tf (t, t + At)} was deduced from it. The 
corresponding {w e (9 y t,t + At)} are shown in figure Fig. 1221 . 

For large angular scales (> 45 degrees), isocontours remain 
open during the whole time range. Large scales have a long coher- 
ence time (~ 5 Gyrs) and are consistent with a 'cosmic tidal field' 
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resulting from the large scale distribution of matter. The latter is 
not expected to evolve significantly with time at our redhsifts and 
the halo's triaxiality is also a fairly constant feature. The innermost 
isocontours are closed around the measurement time ti. Small an- 
gular scales (< 45 degrees) have shorter coherence time (~ 1.5 
Gyr). This is consistent with a contribution to the potential due to 
objects where satellites pass by or dive into the halo and apply a 
tidal field only for a short period. 

This difference between large and small scales can also be 
investigated through the time matrices of the Tf (t,t + At) co- 
efficients (see Fig. 1231 ). The diagonal terms describe the time 
evolution of the angular power spectrum, Tf (t,t) = Cf (t). 
A smooth (resp. peaked) distribution of values around the diago- 
nal indicates a long (resp. short) coherence time. Clearly different 
scales have different characteristic time scales. The non-diagonal 
elements of the quadrupole matrix (£ — 2) decrease slowly with 
the distance to the diagonal while the T20 matrix is almost diag- 
onal. Not surprinsingly, the smaller the angular scale, the smaller 
the coherence time: a small 3D object passing through the sphere 
is likely to have a small angular size on the sphere. 

For a given £ and a given t, the correlations coefficients Tf 
can be fitted by a Lorentzian function defined by 



T e (t,t+At) 



2/tt (At 



W/2) 2 



+<?J e (t).(71) 



where the characteristic time scale, A7V=, is given by gj e (i) and 
the reference time t is equal to qi c (t). Example of fits are shown if 

Fig. 03- 

For example, the time evolution of ATt? = q 2 e (t) is given 
in Fig. I24i for different £ values. Given the error bars, the char- 
acteristic time scales are constant over time (except for the £ = 4 
mode). In the prospect of the regeneration of the potential, the sta- 
tionarity hypothesis can then be considered as valid for most of the 
angular scales. Meanwhile, the I = 4 potential fluctuations display 
a decreasing ATt' with time. The same effect exists at a la level 
for I = 5. One interpretation would be that satellites achieve higher 
velocities along time: for a given typical size a faster satellite would 
spend less time to be accreted and the associated potential would 
be detected on a smaller time scale. This picture is supported by the 
results shown in l5.3.1l where the mean velocity of infalling mate- 
rial increases along time. Another possibility would be that £ = 4 
fluctuations had a longer radial extent in the past. Since there is no 
reasons for potential fluctuations to have such a property, one could 
imagine successive potential fluctuations which overlapped, lead- 
ing to an apparent longer radial extent. This possibility is further 
investigated in the following paragraphs, by comparing the coher- 
ence time variation of the potential fluctuations to the evolution of 
the infall's typical velocity . 
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Figure 22. The angulo-temporal correlation function, w e (8,At) = 
{Si^e-i (CI, t)<5r^,ei (CI + ACl,t + At)). Blue (resp. red) colors stand for 
low (resp. high) values of the correlation. Isocontours are also shown. Large 
angular scale isocontours (8 ~ n/2) have large temporal extent, due to the 
quadrupole dominance over the potential seen on the virial sphere. 



typical time scale AT^e, (w e (AQ.,t, At)} was fitted with a 2D 
function for different values of t. The model used is given by : 



(w e (8,t,At)) = q r(t) + 



q% e {t) 



qj e (t) 



2tt 



(At- q r(t)r+(^: 



x sin(2^(An-gr(t))/g 2 wc (f)) ^ (72) 



AO, 

where the angular dependence is fitted by a cardinal sine function 
while the time dependence is fitted by a Lorentzian function. Ex- 
amples of 2D fits are shown in Fig. I25i . The correlation function 
was also computed using only harmonics with £ ^ 4, 5, 6, 7 and 
the same fitting procedure is applied. The evolution of the resulting 
characteristic time scales AT m i = q± e (t) is shown in Fig. 1261 . 
Bars stand for 3er fitting errors. 

Note that AT w c tends to decrease with time for every trunca- 
tion order. The £ ^ 3 and £ ^ 4 correlation function displays a rise 
of AT,„e before it drops to lower values. Furthermore AT w e tend 
to decreases with £ m in, suggesting that the £ m i n contribution domi- 
nates each w e reconstruction. Correlation functions with £ m i n ^ 5 
show marginal AT w e variation but recall that our time resolution is 
0.53 Gyr, hence any fluctuations on smaller scales should be taken 
in caution. Still, the 0.81 Gyr variation observed for £ ^ 3 between 
t = 5.1 Gyrs and t = 0.8 Gyr is significant and so is the variation 
observed for £ 4 (1.3 Gyr). 



6.1.2.2 Correlations without the dipole and the quadrupole 

In order to focus on the coherence time of small angular scales, the 
correlation function {w e (0, t, At)} was also computed without the 
dipole (|=1) and the quadrupole (£ = 2) component of the poten- 
tial (see also equation 1691 ). The angulo-temporal correlation func- 
tion is shown in Fig. J25I . Again, the isocontours of the correlation 
function are closed around At = 0. This shows that the potential 
on the sphere have a finite coherence time. In contrast to coherence 
times measured on the Tf coefficients, all the angular scales are 
mixed and the typical time scales are those of structures as they 
are 'really' seen from a halocentric point-of view where 'potential 
blobs' appear and disappear on the sphere. To evaluate the related 



6.1.2.3 A longer coherence length As mentioned above, the 
variation of the characteristic time scale can be explained by the 
measured increase in mean velocity. Conversely, the decrease of 
coherence time may be the consequence of smaller potential blobs 
as time passes: a 'large' (three dimensional) potential takes longer 
to disappear than a smaller one. One crude approximation could 
be : 

Li _ V1AT1 
L 2 ~ V 2 AT 2 ' 



(73) 



where L is the radial size of the potential blob, V its radial velocity 
and AT its coherence time. It is assumed that the £ ^ 4 trunca- 
tion is representative of the potential due to infalling objects, i.e. 
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Figure 23. The time matrices of the Tf (t, t+ At) coefficients. Blue (resp. 
red) colors stand for low (resp. high) values of the coefficients. The diagonal 
terms are equal to the angular power spectrum, i.e. Tf (t, t) = Cf (t). 
As can be seen from Fig. 1211 fluctuations observed for T± are within the 
error bars. A smooth (resp. sharp) decrease of Tf (t, t + At) with the 
distance to the diagonal imply a long (resp. short) coherence time. Here, 
coherence time decreases with angular scale. 



AT1/AT2 ~ 1.95. Let us also consider that the radial velocity 
variation is equal to the one measured for the mean velocity of in- 
fall ( see section 15.3.1 \ : V1/V2 = 0.77. Following equation 1731 . 
it suggests that L\ ~ I.5L2, i.e. the radial size decreases with 
time. The same calculation with I ^ 3 leads to L\ ~ I.IL2 : the 
results remain qualitatively the same. In other words, the coher- 
ence length was longer in the past and the velocity variation can- 
not explain the variation of coherence time. The only other way to 
explain a longer coherence length involves potential blobs falling 
successively through the sphere, coming from roughly the same di- 
rection. To induce a decreasing coherence time, these blobs would 
have to be either bigger before or either more numerous. Such a 
crude picture is coherent with the measured decrease of accretion 
with time and the a n isotropic nature of a c cretion by haloes (se e e.g. 
iKnebe et alj l2004 , lAubert et aljEool . lZentner et alJl2005l) ). 
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Figure 24. Top and Middle: the time evolution of the characteristic time 
scale ATt= , obtained by fitting Tf (t,t + At) with equation I7T1 . Sym- 
bols are the measurements while bars stand for 3cr fitting error bars. The 
second order fit of the time evolution of each TJ is also shown. Ex- 
cept for the £ = 4 and (marginally) for the i = 5 modes, no time evo- 
lution is observed. The time resolution is 0.53 Gyr. Bottom: examples of 
Tf (t, t + At) fitted by Lorentzian functions. 
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Figure 25. Top: the angulo-temporal correlation function, w e (8, At) = 
(5 W e] (O, t)<5^e] (fi + Af2,t + At)). The dipole (i = 1) and the 
quadrupole (i = 2) components were removed. Blue (resp. red) colors 
stand for low (resp. high) values of the correlation. Isocontours are also 
shown. The main axes of the 'ellipses' centered on (Af2 = 0, At = 0) 
give indications on the characteristic time and angular scales of ip e (S~l, i). 
Bottom: comparison between the measured 2D correlation function (solid 
lines) and the fit obtained using the equation 1721 (dashed lines). 



6.2 Flux density of mass : w p = pv r 

The mode (C ( T P )) of the distribution of the zc p angular power spec- 
trum is computed using equations 1421 - 1431 . In order to deal with 
adimensional quantities, the reduced harmonic coefficients, 5«, m , 
are defined as: 



at.-. 



4tt( 



('■i. 



UO.Q 



(ooo) (aoo) ' 



(74) 



The accretion contrast, 8[ m ], and the ae tT n coefficients are linked 
by: 



<W°) = {— P ) = 2. 

£ ,m 



Q>£.m 



(75) 
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Figure 26. The cosmic evolution of the potential's coherence time. This 
characteristic time scale is obtained by fitting the 2D correlation function 
w e (8, At) with the function given in equation 1721 . The correlation func- 
tion is computed using harmonics coefficients with I ^ 3 (crosses), I ^ 4 
(squares), I ^ 5 (triangle), I ^ 6 (circle) and i 7 (diamonds). Bars 
stand for 3cr fitting errors. The time resolution is 0.53 Gyr. 



6.2.1 Angular power spectrum 

Given (aoo)(i), the angular power spectrum C^ p {i) is computed 
for each halo. At each time step and for each harmonic order I, 
the C^ p (t) distribution was fitted by a log-normal distribution (see 
Fig.|C4|. The probability distribution of Cj" (t) is weighted as de- 
scribed in appendix 1X1 

The evolution of ((C^ 7p (t)} with time is shown in Fig. 1271 . 
The shape of (C ( T P (t)} remains mostly the same with time and is 
fitted by a simple function : 



<<c7 p >>(*)=«r p (*)+ 



(76) 



(<+ft'(*)) a 

The time evolutions of c7^ p , q^ p and q^ p are shown in Fig. ID2I 
and can be fitted by decreasing exponential: 



Up (t) =/i + fcexp - — 



(77) 



Only the dipole (1=1) harmonic does not fit with the previous 
functional form and is systematically lower than contribution of the 
other harmonics. If particles velocities were measured in an abso- 
lute referential, the dipole strength would reflect the halo's motion 
in the surrounding matter. Also strong mergers may cover a 1 80 de- 
grees angle on the sphere and would contribute to the dipole. Since 
velocities are measured in the halo's rest frame and strong merg- 
ers are excluded, the dipole strength is substantially lowered, as is 
measured. 

The values for h,k,u are given in table ID3l The offset q^ p 
of ((Cj 7p (t)} increases with time. From equation <74t . one can see 
that (Cj 7p (t))) is proportional to the square of the accretion con- 
trast. If the power spectrum experiences a global shift toward higher 
values with time it implies that the accretion contrast increases with 
time. Since the average velocity does not vary strongly with time, 
this suggests that objects are getting denser with time. This effect 
is similar to the global increase of the 3D power spectrum P(k) 
with time due to density growth. Also, the q^ p coefficient is found 
to evolve as q^ p ■ This illustrates the fact that the amplitude of 
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((Cj 7p (i)} remains mainly constant. The q^ p coefficient should be 
seen as a typical scale and varies slightly from q£ p = 6 at z = 1 
to q™ p = 11 at z = 0.1. The {Cj p (t)j becomes marginally 'flat- 
ter' as time passes, implying that small scales contribute more to 
the spatial distribution of uj p (tt, t), consistently with the evolution 
of Ijfiiity' }}. The flat power spectrum measured for zu p on small 
scales suggests that isolated particles contribute significantly and 
increasingly with time. In other words, the accretion becomes low 
enough to be poorly resolved in terms of particles. 



6.2.2 Resolution in mass and particle number 

In order to assess these environments/resolution effects, {C^' p (t)} 
was computed for three different classes of mass (see Fig. [28} at a 
look-back time of 800 Myrs. For the heaviest halos, the power spec- 
trum is peaked toward low £ values. The contribution of large scales 
is quite important. For smaller masses the power spectrum gets flat- 
ter and all scales almost contribute equally for the lightest class of 
mass. Recall that the harmonic decomposition of a Dirac function 
leads to Ce = constant, thus a flat power spectrum indicates that 
isolated particles contribute significantly to the distribution of mat- 
ter on the sphere. The relative behavior of the three ((C J T P (t)} con- 
firms that larger halos still experience important mergers (i.e. on 
large scales) while small ones are in quiet environments at our sim- 
ulation resolution. The effect of the mass resolution on the angular 
structure of accretion was also investigated with two smaller sets 
of simulations : the first one involves 10 simulations with 256 3 par- 
ticles in 50 Mpc/h boxes and the second in 5 simulations with the 
256 3 particles in 20 Mpc/h boxes. The related (C ( T P }} mesured at 
a lookback time of 800 Myrs are shown in Fig. 1291 Here 1532 and 
545 haloes satisfying the conditions described in !4.2l were detected 
in these two additional sets of simulations. For clarity, la error 
bars are shown for the two high resolution measurments while the 
'larger statistics' power spectrum is still represented with 3u bars. 
For large scales (£ < 10), the three power spectra are consistent, 
thus suggesting that convergence was achieved there. On smaller 
scales, the two higher resolution spectra differ significantly from 
the one measured using the other set of simulations (50 Mpc/128 3 
particles) : high £ hold significantly less power. This confirms that 
the lack of resolution tend to overestimates the importance of small 
scales and implies that the study of vj p requires simulations at 
higher resolution in order to understand e.g. the detailed statistics 
of small infalling objects. Interestingly, the two higher resolution 
simulations have identical ((C. p ), given the admittedly large error 
bars. This suggests that statistical convergence at scales £ < 50 
does not requires extremely resolved simulations and simulation 
boxes with a mass resolution only 8-10 times greater than those 
used in this paper should suffice. 

Finally, it clearly appears from Fig. 1271 and Fig. <D2t that the 
angulo-temporal correlation function related to {C e p }(t) are dom- 
inated by the overall shift of the angular power spectrum towards 
higher values and consequently no coherence time should be de- 
tectable. It is shown in appendix|E|how the evolution of {C e p }(t) 
is related to mass biases and possible resolution effects. The pre- 
vious measurements on the potential were not sensitive to these 
effects because of the smoother nature of the field. 

In appendix |E| the measured secular evolution is avoided by 
using an alternative definition of {C^ p }(t). It is found that the 
angular power spectrum of zu p can be fitted by a simple power-law 
(see Fig. lEH at every time : 
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Figure 27. The average angular power spectrum, (Cj 7p )(t), at t = 
0.8,3.5,5.7,7.8 Gyr (symbols). {Cj p ){t) is taken as the mode of the 
log-normal function used to fit the Ci distribution. Bars stand for 3<r errors. 
For a given I the corresponding angular scale is ir/£. (C e p )(t) may be 
fitted by a generic model given by equation 1761 (solid line). 
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Figure 28. The average angular power spectrum, (C^ p )(t), at t 
Gyr (symbols) for three different class of masses. (CY P ){t) is taken as 
the mode of the log-normal function used to fit the C^ p distribution. Bars 
stand for 3<r errors. Mass are expressed in solar masses. For a given I the 
corresponding angular scale is n/£. The three measurements are fitted by 
equation 1761 (solid line). The power spectrum gets natter for small halos. 
Accretion by small halos is dominated by small objects or even isolated 
particles. 



The corresponding angulo-temporal correlation function is given in 
Fig. lE3l As expected, there is a shorter coherence time for the vj p 
field than for the potential, because of the 'sharper' nature of the 
former. Overall, these results strongly suggest that class of mass 
and resolution biases should be systematically investigated beyond 
what was shown here. 
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Figure 29. The average angular power spectrum, {C ( p }{t), at t = 0.8 
Gyr for three different mass resolution in the simulations. {C™ p )(t) is 
taken as the mode of the log-normal function used to fit the C^ p distri- 
bution. For a given i the corresponding angular scale is n/£. Circles stand 
for the measurements performed of the main set of simulations (50 Mpc/h, 
128 3 particles) while error bars stand for 3 — a errors. Square and tri- 
angles stand for measurments performed on simulation with higher mass 
resolution, resp. 50 Mpc/h-256 3 particles (1532 haloes analysed) and 20 
Mpc/h-256 3 particles (545 haloes analysed). In these two cases, error bars 
stand for 1<7 errors. 



7 A SCENARIO FOR THE ACCRETION OF A TYPICAL 
HALO AT Z 1 

Let us draw a summary of the previous results, in order to get 
a synthetic picture of the flux properties at the virial radius. A 
typical halo in our sample has a mass of 1O 13 M0 and a radius 
^200 ~ 500kpc at z=0. It is embedded in a quasi-stationary gravi- 
tational potential, ip e . Such a potential is highly quadrupolar and is 
likely to be induced by the large scale distribution of matter around 
the halo. The halo accretes material between z = 1 and z = at 
a rate which declines with time. At high redshift, only accretion is 
detected at i?200 • It is mainly radial and occurs at a velocity close 
to 75% times the circular velocity. As time advances, accretion of 
new material decreases, while outflows become significant. Out- 
flows occur at lower velocities and on more circular orbits. A frac- 
tion of the outflowing component is due to a backsplash population 
made of material which already passed through the virial sphere. 
Another fraction of outflows corresponds to 'virialized' material of 
the halo which goes further than R 2 oo and is being 'cancelled' by 
its infalling counterpart. The clustering on the sphere of the gravita- 
tional potential drifts toward smaller scales, while the clustering of 
matter follows marginally the same trend at our level of resolution. 
It reflects the increasing contribution with time of weak/diffuse ac- 
cretion, poorly sampled at our resolution. In parallel, the coherence 
time of potential fluctuations is found to be decreasing with time by 
the halocentric observer. This decrease may be related with the ac- 
cretion of satellites, where objects were numerous enough to 'over- 
lap' in the past, which implies that accretion occurs mainly in the 
same direction on the virial sphere. 

This scenario seems consistent with most of the past studies 
made on the subject using the full 3D information contained in 
simulations. The decline of accret ion rate has been already mea- 
sured by e.g van denBoschir2002bl) even though our measurements 



are in a slight quantitative disagreement. The 'rebound' of mat- 
ter through the Viri al sphere has alread y been measured by e.g. 
iMamon et ail 120041) or lGill et all <2004 . Furthermore, the veloc- 
ity bimodality is recovered together with the circularisation of or- 
bits measured bv lGill et alJ (2004) in high resolution simulations. 
Finally, the variation of the potential's coherence time is found to be 
related to the aniso t ropy of accretion, al rea dy demonstrated by e .g. 
iKnebe et alJ <2004 . lAubert et allJiooH) or lZentner et"al]l2005l) . 



8 SUMMARY & DISCUSSION 

This paper, the second in the series of three, presents measurements 
of the detailed statistical properties of dark matter flows on small 
scales (^ 500 kpc) in the near environment of haloes using a large 
set of ACDM cosmological simulations. The purpose of this in- 
vestigation was twofold: (i) characterize statistically (via one and 
two points statistics) the detailed (angular and kinematic) incom- 
ing fluxes of dark matter entering the virial sphere of a biased (de- 
scribed in Section l4~2l sample of halos undergoing minor mergers 
for the broader interest of astronomers concerned with the environ- 
ment of galaxies; (ii) compute the first two moments of the linear 
coefficients, c n , (resp. 6 n ) of the source term (resp. external poten- 
tial) entering equations 1 1 3i - i 141 (paper I). 

We concentrated on flows at the haloes' virial radius while de- 
scribing the infalling matter via fluxes densities through a spherical 
shell. In parallel, we measured the statistical properties of the tidal 
potential reprojected back onto the boundary. The statistical one 
and two point expectations of the inflow were tabulated both kine- 
matically and angularly on the -R200 virtual sphere. All measure- 
ments were carried for 15000 haloes undergoing minor (as defined) 
mergers between redshift z=l and z=0. The two points correlations 
are carried both angularly and temporally for the flux densities and 
the tidal field. We also provided a method to re-generate realization 
of the field, via Equation I38i and <F6t ). 

We briefly demonstrated how a perturbative description of the 
dynamics of halos can propagate the statistical properties of envi- 
ronments down to the statistical properties of the halo's response. 
The description of the environment involved the projection of the 
potential and the source on a basis of functions. This basis allowed 
us to decouple the time evolution from the angular and velocity 
dependance of these two quantities. Hence, the accretion and the 
tidal potential were completely described by the projection coeffi- 
cients and their statistical properties, which depend on time only. 
We also discussed how the flux densities of matter, momentum and 
energy can be related to the source and its expansion coefficients. 
We restricted ourselves to the one and two statistical description of 
the tidal fiel d ?/> e and the flux density of mass vo p and postponed 
to paper III jAubert & Pichonl feood) ). the full description of the 
higher moments. Since these measurements will be used as an en- 
try to a perturbative description of the inner dynamics of haloes, 
only objects with quiet accretion history were selected as discussed 
in Section l4~2l Throughout this biased sample of haloes we made 
statistical measurements of the kinematic properties of accretion 
and derived results on the following quantities : 

• the evolution of the accretion rate at the virial sphere: the net 
accretion is found to decrease with time, probing both the increas- 
ing contribution of outflows and the decline of strong interactions. 

• the evolution of the net velocity distribution of the accretion : 
infall exhibits a typical velocity of 0.75V C . A backsplash compo- 
nent is detected at recent times with a significant outflowing com- 
ponent at a lower velocity (~ 0.6V C ). 
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• the evolution of the impact parameters/incidence angle distri- 
bution of the infall. The infall is found to be mainly radial while 
outflows are on more circular orbits. 

• the angulo-temporal two-point correlation of the external po- 
tential on Virial sphere. The potential appear to be mainly domi- 
nated by a strong and constant quadrupole. The coherence time of 
smaller angular scales provide hints of an anisotropic accretion. 

• the angular power spectrum of accreted matter. The cluster- 
ing is dominated by small angular scales, possibly at the resolution 
limit. 

• the angulo-temporal correlation of the flux density of mass. 
The coherence time appears shorter than for potential fluctuations, 
as expected. 

These results can be interpreted in terms of properties of accreted 
objects or of smooth accretion and are coherent with previous 
studie s (e.g.lGhigna et alj Jl998h.lMamon et alj j2004h.lGill et all 
2004. Ivan den Boschl <2002bh. iKnebe et all <2004. lAubert et all 



2004lBensoiI (2005)). These studies were mainly focussed ont 



the properties of accreted satellites. Properties of accreted subha- 
los could also be directly derived from these results, once a clear 
definition allow us to distinguish structures within the general flow 
of matter. Substructures are expected to form a distinct "phase" of 
the accreted fluid: for instance the velocity dispersion is expected 
to be quite different in compact objects than in the smooth accret- 
ing component. This phase separation will be assessed in paper III 
where a systematic comparison of the current approach with an 
analysis in terms of pre-ide ntified satellites will be carried, in the 
spirit o flAubert et al.li2b04l) . The contribution of outflows, the lack 
of a standard definition for subhaloes, resolution issues and the fact 
that properties are measured at one radius (which could be statis- 
tically propagated toward inner regions) are all issues which must 
be assessed before a complete and rigorous comparison can be per- 
formed. Yet, the current agreement between a fluid description of 
the environment and these above mentioned published results are 
clearly encouraging hits for the reliability of the method presented 
here. 

These kinematic signatures provide insights on the processes 
which occur in the inner regions of haloes. In particular, the kine- 
matic discrepancies between the different components of the mass 
flux should be understood in terms of dynamical friction, tidal 
stripping or even satellite-satellite interactions within the halo. The 
kinematical properties of accreted matter may be transposed to 
the kinematical properties of satellites observed around galaxies. 
Newly accreted material exhibits kinematic signatures (radial, high 
velocity trajectories) different from the ones measured for matter 
which already interacted with the halo (tangential, low-velocity or- 
bits). Admittedly, it is not straightforward to apply directly these 
results to the luminous component (see paper I for a discussion 
of thresholding), and to see how projection effects may affect the 
distributions. Nevertheless, the corresponding observationnal mea- 
surements on satellites should provide information on the past his- 
tory of these objects. 

As discussed in paper I, these measurements can be used as 
an entry to the perturbative theory of the response of the open halo. 
Phenomena related to accretion can be consistently assessed via 
this framework: dynamical friction, tidal stripping and phase mix- 
ing. With the statistical description of the tidal field presented in 
this paper only, we may already implement the theory presented in 
paper I in the regime of pure tidal excitation. The complete knowl- 
edge of the source (which will be completed in paper III) should 
considerably extends the realm of application provided by this the- 



ory. Specifically, we have shown in paper I that the internal dynam- 
ics of sub structures within galactic haloes (distortion, clumps as 
traced by Xray emissivity, weak lensing, dark matter annihilation, 
tidal streams ..), and the implication for the disk (spiral structure, 
warp, etc..) could be predicted within this framework. Conversely 
the knowledge of the observed properties of a statistical sample of 
galactic halos could be used to (i) constrain observationnaly the sta- 
tistical nature of the infall (ii) predict the observed distribution and 
correlations of upcoming surveys, (iii) time reverse the observed 
distribution of clumps, and finally (iv) weight the relative impor- 
tance of the intrinsic (via the unperturbed distribution function) and 
external (tidal and/ or infall) influence of the environment in deter- 
mining the fate of galaxies. 

The current measurements reduce the degree of freedom that 
still exists in the setting of numerical experiments in a galactic con- 
text. For instance, given that the structure of the external tidal field 
is found to be simple, it can easily be modelled as an external com- 
ponent in numerical simulations (or even in analytical studies). It 
would provide a simple but statistically relevant contribution of the 
large scale structures to the dynamical states of haloes.The tem- 
poral coherence of the first £ > 2 angular harmonics of the tidal 
field should allow one to draw more accurate representations of ex- 
ternal contribution to the field that would include the fluctuations 
due to smaller structures. The kinematics of accretion is not ran- 
dom as well and the distribution of velocities at R200 follows a 
gaussian-shaped curve which caracteristics evolve with time and 
exhibit a certain distribution of the impact parameter. These re- 
sults put prescriptions that could be used to generate encounters 
between satellites and galactic disks that follow the ones measured 
in cosmological simulations at large radii. We also presented first 
constrains on the angulo-temporal correlation function of the ac- 
cretion. Even though it is not completely clear yet how resolution 
will eventually affect these results, such functions contains some 
glimpse of informations regarding the angular distribution of en- 
counters with external systems but also regarding the frequency of 
accretion events. This frequency can also be probed by the tem- 
poral coherence of the fluctuations in the tidal field. The apparent 
contradiction that exists between the observed number of discs and 
the predicted large number of mergers may be solved by a better 
knowledge of the frequency of the latter: it may be low enough to 
solve this contradiction. In this context, simulations of successive 
mergers between a galaxy and satellites should be consistent with 
large scale simulations an we provide first constrains on the rate of 
minor encounters at the halo's outer boundary. 

As argued in lPichon & AuberJ l200(j|) . we emphasize that an 
a priori discrimination between 'objects' and diffuse matter may 
not constitute the best way to describe accretion: it is not clear that 
luminous matter is always attached to dark matter overdensities, 
there is no unambiguous definition of substructures and their state 
change with time under the influence of tidal shocks or dynamical 
friction. The generation of objects that follow the current results is 
admittedly the easiest way to proceed but should be followed by 
a more general description of matter in terms of 'fluid approach', 
where 'objects' only constitute a specific phase of such a fluid. The 
statistical measurements on both ip e and tjj p allow the regeneration 
of synthetic environments. Knowing the average evolution and the 
angular power spectra of these quantities, the generation of spher- 
ical maps in the gaussian regime is straightforward. We describe 
such a regeneration procedure in appendix|F| Such maps would ef- 
ficiently provide realistic environments of halos, consistent (up to 
two-points statistics) with those measured in cosmological simula- 
tions and could be 'embedded' into simulation of galaxies. Again, 
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virialised structures would be naturally included (since they have 
their own statistical signature on the virial sphere) without relying 
on any adhoc prescription on their nature. 

E xtensions to non gaussian fie lds are also possible (follow- 
ing e.e. lContaldi & Magueiiol feoOlh ) but would rely on higher or- 
der correlations. It was assumed throughout this investigations that 
fields could be approximated as gaussian fields, fully described by 
their two-points statistics. Yet a simple visual inspection of w p 
maps reveals that they are not strictly gaussian, a finding confirmed 
by preliminary analysis of their bispectrum. Furthermore, paper I 
demonstrated that a dynamical description which takes into account 
non-linear effects, such as dynamical friction, requires higher-order 
correlations. Therefore, extensions to non-gaussian fields are in or- 
der in the long run. 

It should again be emphasized that some aspects of the present 
work are exploratory only, in that the resolution achieved (Mhalo > 
5 ■ 1O 12 M0) is somewhat high for galaxies. In fact, it will be 
interesting to confirm that the properties of infall do not asymptote 
for lower mass (Mh a i < 5 • 1O 12 M0) together with the intrin- 
sic properties of galaxies. In addition a systematic study of biases 
induced by our estimators of angular correlations should be con- 
ducted. For a fixed halo mass, our lack of resolution implies that 
we over estimate the clumsiness of the infall. 

As demonstrated in Section l6,2,2l the limited resolution (both 
spatially and in mass) of our simulations appeared to be an issue 
for some of the results presented here (e.g. the angular power spec- 
trum of vj p ). Systematic use of higher resolution simulations (in the 
spirit of section 6) will be required to fully assess these limitations. 
In particular, a fraction of the accretion detected as a weak/diffuse 
component may be associated to unresolved objects; the influence 
of small mass satellites should therefore be explored. With the 
prospect of deducing the properties of galaxy from halos environ- 
ments, lower mass haloes are more likely to host only one galaxy, 
making them more suitable for such a study. Cosmological simu- 
lation of small volumes also tend to prevent the formation of rare 
events which may be relevant for the representativity of the study: 
for example, some disks seem to indicate that they were formed 
in 'very quiet' environments. The right balance between resolution 
and volume should be found. Aside from these biases induced by 
simulations, we also introduced selection criteria on both the mass 
or the accretion history of haloes and the influence of these arbitrary 
choices on our statistical distribution should be assessed precisely. 

Eventually using hydrodynamical codes which include bary- 
onic effects in simulations and introducing the physics of gas in our 
model, we would construct a complete semi-analytic tool to study 
the detailed inner dynamics of galaxies. 
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In order to set the maximal order of this expansion, ^ max , we 
compared the two-point correlation function of the spherical field 
pv r (Cl) to the one inferred from harmonic coefficients, a £m , and 
power spectrum Ci (see equations (equation <42t ) and (equa- 
tion <43H ), Within a set of randomly distributed particles, let 
d P oisson(#) be the probability of finding two particles with an angu- 
lar separation 9. If d(9) is the same probilty for a given distribution 
of particles then its two points correlation function is defined 



d(9) 



dm 



(Al) 



(0) • 

The co rrelation functio n and the ae m coefficients are related 
by (e.g. lPeacocklJ 1999h ) : 



£(6>) = Y^C e (2£+l)P e (cos( 



(A2) 



Figure Al. The average angular two-point correlation function, 
of the advected mass spherical field pi> r (f2). The correlation function is 
shown as a function of the angular distance on the sphere 8 given in radi- 
ans. Using the positions of accreted particles around 200 halos at z=l, the 
average correlation function can be computed (dots). Lines represent the 
correlation function deduced from the harmonic coefficients of the pv r (fl) 
fields around the same 200 halos, with lmax= 10, 30, 50, 60. The conver- 
gence is ensured for lmax ^ 50. 
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APPENDIX A: HARMONIC CONVERGENCE 

As explained in section |3~T1 the angular dependence of the source 
function is expanded over a basis of spherical harmonics Yi m (fl). 



where 9 is an angular distance on the sphere and Pi(x) a Legen- 
dre function. The average angular correlation function 
is defined as: 



■^2n 2 p £ p (9), 



(A3) 



where (, p (9) is the two-point correlation function of the p-th halo 
computed using n p particles passing through the virial sphere. 
From a set of 200 halos extracted from a simulation, we computed 
{£,{9)} using different values for i! m ax (see Fig. lAU . From the same 
set of halos, we also computed the average two point correlation 
function directly from the particles positions using equation <AU . 
From Fig. lAU it clearly appears that ((£(#)}} pair has not con- 



verged for 



< 30. For . 



50 the actual two point cor- 



relation function is well reproduced. Since no real difference can 
be distinguished between £ max = 50 and i! m ax = 60, we chose to 
limit the harmonic expansion of the source term to £ ^ 50. Note 
that the truncation in ^ ma x defines an effective resolution beyond 
which the distribution is effectively coarsed grained. 



APPENDIX B: ANGULO-TEMPORAL CORRELATION 

Let us consider a spherical field X(tt, t) which can be expanded 
over the spherical harmonic basis: 



(in 



(Bl) 



The correlation w x between two successive realizations of X is 
defined as: 



W X (Sl,Cl',t,t') EE (X(il,t)X(Sl',t')) 



(B2) 



where (■) stands for the statistical average. If X(tt, t) is isotropic, 
the correlation should not depend on 17 or fi' but only on the dis- 
tance 9 between the two points. It implies that w can be expanded 
on a basis of Legendre Polynomials Pl{u) '■ 

w x {ft,ft',t,t') = w x (9,t,t') = ^(2L+l)T L P L (cos((9)).(B3) 

L 

How are Tl and xi m related ? Rewriting equation <B2> as: 

w x (9,t,t') =^^(x« m (i) a: l, m ,(t , )>^m(n)Y; m ,(n') J (B4) 

£m I'm' 

one can write: 
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J <mdn'Y; imi (n)Y t2m2 (n')w x = (x <imi (t)ij i , m3 (t')).(B5) 

Meanwhile, assuming isotropy, one can also write 

J <mdn'Y; imi (n)Y e2m2 (n')w x = 

^(2L + l)Ti f dSldfl'Y t * imi (U)Y t2rn2 (U')P L (co S (e)) 

L J 

= ]T(47r)T L J dndn%\ mi {n)Y l2m2 {n')Y LM {n)Y* LM {n') 

LM J 

1 SLt 2 SMm 1 SMm 2 i (B6) 

LAI 

where Pl is expressed in terms of spherical harmonics using the 
spherical harmonics addition theorem. In the end, we get: 



10* 6 10* 6 
0. 2. 4. 6. 8.0. 2. 4. 6. 




300 



T e = -r-{xt m (t)xi m (t')). 

For a given realization of X(fl, t), Te can be estimated by: 
1 1 



4tt 21 + 



(B7) 



(B8) 



APPENDIX C: DISTRIBUTIONS 

In this appendix, we present the various distributions fitted either by 
normal or log-normal PDF. For each quantity, the mode (or most 
probable value) of the distribution has been obtained from these 
fits. The gaussian distribution is defined by: 



N(x) 



A 



exp 



(s ~ A*) 5 

2cr 2 



(CI) 



while the mode is equivalent to the mean fx. The log-normal distri- 
bution is given by 



LN(s) = 



.4 



exp 



(logQ//*)) 5 
2cr 2 



(C2) 



while the mode is given by /iexp(— a 2 ). The different fits men- 
tioned in the main text are described in the following figures: 

• Fig. IC II shows the distributions of the harmonic coeffi- 
cient boo (i) which is proportional to the potential averaged on the 
sphere. It is expressed in units of GM/R where M is expressed in 
1O 1O M , R in kpc/iT 1 and G=43007 in internal units. 

• Fig. <C2> shows the distributions of the external potential's 
power spectrum for four different harmonic £ = 2, 5, 10, 20 at 
t=1.3 Gyr. The distribution has been fitted by a log normal func- 
tion. 

• Fig. <C3t shows the distributions of the mean flux $ M (t). The 
mean flux is proportional to the harmonic coefficient aoo- The dis- 
tribution is fitted by a normal distribution. The normal model agrees 
well with the measured distribution at recent times but fail to re- 
produce the outliers tail at high redshift. Consequently the mode 
position is underestimated at these times. 

• Fig. <C6t shows the distributions of one of the coefficient in- 
volved in the computation of the velocity distribution 4>(v). Four 
different times are being represented. The coefficient distribution 
has been fitted by a gaussian distribution. 

• Fig. <C4> shows the distributions of the power spectrum val- 
ues CT P for four different harmonic order I. The fits were made 



^Jll Mil 111] Mil II II 1 1 1 1 Mil |HU£^ 

E ffi t=7.8 Gyri| 


1 1 111 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 [L^ 

: A t=5.6 Gyr % 


\k t=2.9Gyr -if 

m\\ nil lirnftrirhfrrh rti mm imf 1 


\ t=.27 Gyr] 

tiiLmi^rit^^ 



0. 2. 4. 6. 8.0. 2. 4. 6. 8. 

10* 6 10* 6 

-boo 

Figure CI. The probability distribution of boo at four different times. This 
coefficient is proportional to the external potential averaged on the sphere. 
The log-normal fit is also shown. 
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Figure C2. The probabiUty distribution of Cf° the for £ = 2,5, 10, 20 at 
the look-back time t = 1.9 Gyr. Note that the X-axis is sampled logarith- 
mically; the corresponding log normal fit of the mode is also shown. 



at t = 1.9 Gyr. This distribution has been fitted by a lognormal 
distribution. 

• Fig. IC5I shows the distributions of the power spectrum values 

C, p for four different harmonic order I. The fits were made at 
t = 2.95 Gyr. This distribution has been fitted by a lognormal 
distribution. 



APPENDIX D: FITS AND TABLES 

In the main text, some statistics are fitted by simple laws and the 
time evolution of these statistics can be described by the time evo- 
lution of the fitting parameters. In this appendix, we present the 
time evolution of these parameters: 
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Figure C3. The probability distribution of the mean flux <£ M (i) for four 
different times. The normal fit is also shown. 



Figure C5. The probability distribution of C e " the for I = 3, 10, 25, 40 
at the look-back time t = 1.9 Gyr. Note that the X-axis is sampled loga- 
rithmically; the corresponding log normal fit of the mode is also shown. 
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Figure C4. The probability distribution of the Cg for i = 3, 10, 25, 40 at 
the look-back time t = 1.9 Gyr. Note that the X-axis is sampled logarith- 
mically; the corresponding log normal fit of the mode is also shown. 
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Figure C6. The probability distribution of the c^o ' or me l°°k-back 
times t = 1.9, 3.5, 4.5 and 5.7 Gyrs. Note that the X-axis is sampled lin- 
early; the corresponding fit of the mode is also shown. 



• Fig. JD1> and Table lD ll presents the time evolution of the ve- 
locity distribution of accretion. It can be fitted by two gaussians 
drifting toward higher velocities with time. 

• Table lD2l summarize the time evolution of the distribution of 
incidence angles. This distribution can be fitted by a gaussian (for 
the inflows) and a linear relation (for the outflows). The fitting pa- 
rameters present a linear evolution with time. 

• Fig. jD2j and Table IMl summarize the time evolution of the 
angular power spectrum of zu p . The power spectrum can be fitted 
by equation iJ6\ where fitting parameters decreases exponentially 
with time. 



</l.l 


-0.028 


0.912 


9i,2 


0.0059 


0.084 


<?i,3 


2.1 10~ 8 


9.03 10~ s 


Qo,l 


-0.026 


0.54 


1o,2 


-0.028 


0.23 


(*) 
<3 


1.79 10~ 8 


-1.22 10~ 7 



Table Dl. Fitting parameters for the time evolution of Gaussian coeffi- 
cients, q(t) = m X t(Gyr) + n. These coefficients, together with equa- 
tion 1611 allow us to compute the time evolution of (<p(v, t)). (*) By defi- 
nition, q Di 3 ^ 0. 
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Figure Dl. time evolution of the fitted coefficients (given in table IdTT 
of the PDF {tp(v, t)) for the parameterization suggested in equation Id 1 1 
(symbols). Bars stand for the 3a errors on the coefficients determination. 
Qi k>k = 1) 2,3 are respectively the mean, the r.m.s and the amplitude 
of the positive Gaussian (i.e. infalling component). q j., k = 1, 2, 3 are 
the corresponding values for the negative Gaussian (i.e. outflowing compo- 
nent). The trend is accurately fitted by a linear relation (red dashed line). 
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Table D2. Fitting parameters for the time evolution Pfe(t), k = 1, 2, 3 pa- 
rameters. p(t) = m X f-(Gyr) + n. Together with equation 1651 . these 
coefficients allow to predict the time evolution of (i?(ri , i)). 
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Figure D2. Time evolution of the coefficients for the (C^ p )(t) model 
(see equation il(>Y )(svmbols). The three time evolutions may be described 
accurately by decreasing exponentials (see equation 1771 ) (solid lines). 
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3.153 
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VJ p 


5.4105 


9.44 
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Table D3. Fitting parameters for the time evolution of (Ci} TUp model's 

coefficients (see (equation 1761 )1 q^ p (t) = h + kc with lookback 
time t expressed in Gyr. 
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Figure El. The angular power spectrum of the external potential, 
(C e rUp )(t), at four different lookback times (symbols). Harmonics coef- 
ficients were normalized using equation IE2I . halo by halo. {C^' p }{t) is 

taken as the mode of the log-normal function used to fit the C ' ^ p distribu- 
tion. Bars stand for 3cr errors. For a given I the corresponding angular scale 
is it /I. The power spectra maybe fitted by a generic model given by equa- 
tion (E3j (solid line). Unlike (C™")(t), (C^ p )(t) remains the same 
with time because of a different normalization. 

APPENDIX E: ALTERNATIVE CONTRAST AND 
ANGULO-TEMPORAL CORRELATION 

It appears clearly from the time evolution of (Cj 7p }} (t ) presented in 
the main text, that no angulo-temporal correlation can be computed, 
since it would be completely dominated by the secular evolution of 
the power spectrum. For this reason, an alternative definition of the 
flux density contrast has been used: 



sum 



or in terms of harmonic coefficients: 



4tt( 



Sio). 



(El) 



(E2) 



Using this new definition, the fluctuations of vj p on the sphere are 
measured relative to the mass flux of each individual halo. The main 
drawback of such a definition is that it couples the etoo and the at m 
coefficients from the beginning. For example, the resimulation of 
such fields would require the knowledge of conditional probabili- 
ties, i.e. what is the distribution of ai m for a given value of aoo, 
while the previous definition (given by equation <74t ) only acted as 
specific choice of (non constant) units. Furthermore this new defi- 
nition is more sensitive to outliers with low aoo values. 
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Still, the angular power spectrum {C/^''}(t) of 5'^ j(O) is 
much more regular than the one obtained from the previous defi- 
nition. Its overall amplitude remains constant over the last 8 Gyrs, 
while its shape seems to be less dominated by small scale contri- 
butions. This alternative power spectrum is well fitted by a single 
power law : 

{C' t w >)(t)=0.75l- 1M , (E3) 

for the whole time range covered by the current measurements. 
This constant shape suggest that harmonic coefficients scale like 
cioo, i- e - the mass flux. Such a scaling is not obvious, since a 
strongly clustered zu p field may coexist with a nil net flux (i.e. 
aoo ~ 0). It also implies that the evolution measured on the previ- 
ous definition of the power spectrum, {C e p )} (t), is more related to 
the evolution of the average flux (traced by aoo) than to the mod- 
ification of the fluctuations amplitude (traced by the others ai m ). 
Still, the evolution of (aoo) spans over one magnitude while the 
evolution of {C™ p }(t) spans over several order of magnitudes: 
this strongly suggests that two different populations of haloes con- 
tribute to the two type of power spectrum. In Fig. <E2> . the scatter 
plot of Cao and C 40 as a function of aoo shows that the haloes 
which experience strong accretion dominate the peak of the C40 
distribution while haloes with low accretion dominate the peak of 
the C40. Furthermore C40 does not scale anymore like aoo as it 
drops below some level, providing hints of resolution and isolated 
particles effects. To conclude, 5i w 1 (SI) appears as better way to 
rescale the fluctuations' amplitude since it provides a more regular 
behavior of the power spectrum, but it is a more complex quantity 
to manipulate. Meanwhile, <5r ro 1 (SI) is probably the correct way to 
to proceed but is clearly more sensitive to resolution effects, which 
should be assessed with bigger simulations in the future. 

Since the behavior of 5[ ro ] (SI) is more regular than the previ- 
ous contrast definition , the angulo-temporal correlation function of 
the flux density of mass has been computed from this new defini- 
tion. Since this definition is different from that used for the poten- 
tial, we restrict ourselves to a qualitative description. The correla- 
tions are given in Fig. IE3i . Clearly, the correlation is more peaked 
around At = and more generally w^ p is sharper than to* 1 . Note 
that no multipole i has been removed during the computation of 
w mp , implying that the quadrupole effect measured in the poten- 
tial correlation is simply not detected for this field. This strongly 
suggest a large-scale 'cosmic' origin for the quadrupolar tidal field 
rather than an artifact of the spherical intersection of an ellipse. Fur- 
thermore, the correlation time is smaller than the one measured for 
the potential, even compared to the correlation time of the potential 
without the 1 = 2 component. This is coherent with the fact that 
density blobs should be 'sharper' than potential blobs as they pass 
through the sphere. 



APPENDIX F: RE-GENERATING GALACTIC TIME 
LINES 

Given the measurements given in Section [5] and Section |6| let us 
describe here how to regenerate realizations of the history of the 
environment of haloes, first for the tidal field only, and then for the 
full accretion history. 
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Figure E2. Scatter plots of the power spectra C^f (Top) and C'^ p (Bot- 
tom) as a function of <2q . The four colors stand for different lookback 
times: t= 7.8 Gyrs (red), 5.6 Gyrs (green), 4. Gyrs (yellow) and 2.8 Gyrs 
(blue). The monopole aoo scales like the integrated flux of matter. The 
quantities C 40 P and C 40 p differ by the normalization applied to the har- 
monic coefficients a^ m . See the main text for more details. 



Fl Re-generating tidal fields 

Let us first focus on the generation of the potential tidal field gener- 
ated by fly by's, hence neglecting the influence of the infall through 
the virial sphere. 5 

First we consider two time variables: T (the 'slow time') 
which describes the secular evolution of the field and tf (the 'fast 
time') which describes the temporal evolution around a given value 
of T, hence describing high frequencies variations. We assume that 
correlations exist only on small time scales, while variations on 
the 'slow time' scale describe secular drifts. Therefore, the field's 
regeneration should include both correlations on short time scales 
and long term evolution. 

Let us call = ip m]UJ (T) the temporal (with respect to the 



5 note that this assumption is not coherent with the way the measurements 
are carried, since it implies that the infalling material somehow disappear 
after crossing -R200 ■ 



32 Dominique Aubert and Christophe Pichon 




angular scale (rad) 

Figure E3. The angulo-temporal correlation function, w ^(6, At) = 
{<Sr OT i (fi, t)5[ ro ] (fJ + AO, t + At)). Blue (resp. red) colors stand for low 
(resp. high) values of the correlation. Isocontours are also shown. Large an- 
gular scale isocontours (Af2 ~ 7r/2) have large temporal extent, due to the 
quadrupole dominance over the potential seen on the virial sphere. 



fast time) and angular Fourier transform of the potential, at fixed 
slow time, T. The probability distribution of "P, is given by 



exp 



■C7 1 - 



(*>) 



\l/2 



det 



p(*) 

(2tt) 

where the variance reads 

C*(T) = ((* m , w - <* m , w )).(* m , u - (*m, W ») 
and the mean field obeys 
(*>(T) = (* m , w ) . 



(Fl) 



(F2) 



(F3) 



Since the potential is isotropic, ($?)(T) is essentially zero (see 
also Fig. IFH . while stands for the angular power spectrum 
described in section lo"Tl Let us call {"I'm, u(T;)}i^ at, the set of 
sampled 4/ a fixed slow time, T 4 . Relying on a linear interpolation 
between two such realizations, the corresponding external potential 
reads in real space: 



il> e (t, fl) = ^ ^ ^""^ J dui exp (lm ■ f2 + luit) /C™ (t, cu) 

i m 

where the kernel, JCf 1 (t, uj), reads 



(F4) 



icT{t,u) = 



,{T 



t - Ti 



+i) 



Ti 



Ti+i — Ti 

Q(T i+ i - t)Q{t - Ti) , (F5) 

recalling that Q(x) is the Heaviside function. 

To sum up, the computation of ^m,^ following IF1I and IF3I 
ensure that correlations on short periods are reproduced, while the 
interpolation procedure allows to take in account the long period 
evolution of the field. This procedure can be repeated for an arbi- 
trary number of virtual tidal histories. 

F2 Re-generating tidal fields and infall history 

Let us assume briefly that the fields are stationary both in time and 
angle, and that their statistics is Gaussian. As shown in Fig. lFll this 



assumption is essentially valid for the expansions of the potential 
and the flux density of mass. Let us call the eleven dimensional 
vector II(t) = (vj p , zu pv , zu prTtT ,tp e ); and Ilm,^ the temporal and 
angular Fourier transform of the fields. The joint probability of the 
field, p(n m ,w), reads 



p(n m ,w) = 
where 



exp 



i(n-<n» .c^.(h-(h)) 



Cn 

and 



(n: 



(2*) 

,«» ■ (n 



11/2 det 1/2 |C, 



m ,u 



<nU)))] 



<n> = [(nL 



(F6) 



(F7) 



(F8) 



Since these fields are mostly isotropic, their expansions coefficients 
are nil on average. Hence, the quantity Cn stands for the angular 
power spectrum of the eleven fields. For respectively the potential 
and the flux density of mass, its temporal evolution is described in 
section lo"Tl and l6,2'l . These measured power spectra are sufficient to 
generate environments restricted to the flux density of mass and the 
potential. We emphasize that these two fields would not be coherent 
if no-cross correlations is specified. These cross-correlations and 
the power spectra for higher moments of the source are postponed 
to paper III. 

Assuming the full knowledge of these 11 fields and their 
cross-correlation, it is therefore straightforward to generate for each 
(m, w) a eleven dimensional vector which satisfies equation IF6i . 
and repeat the draw for all modes (both lj and m). Inverse Fourier 
transform yields II(t). Once 11(f) is known, we can regenerate the 
whole five dimensional phase-space source as a function of time 
via equation I38i . This process can also be repeated for an arbitrary 
number of virtual halo histories. The assumption of stationarity in 
time can be lifted following the same route as that sketched in Sec- 
tion lFll tsee Equation IF5I ). 



APPENDIX G: FROM EXPANSION COEFFICIENTS TO 
FLUX DENSITIES 

Gl From expansion coefficients to advected momentum 

The phase space distribution of advected momentum is given by: 

uj pv (fi,v,r,t) = s e (ci,v,r,t)v. (Gi) 

= ™ m >(t)g a (v)Y m (n)Y m ,(r)v,(G2) 

a ,m,m' 

where the velocity vector may be written as a function of spherical 
harmonics: 

v = -v^i-Y^W + Y^TVee 



v^-YiKDer. 



Then, one can write 
uj pv (Cl,t) = 



dvdrv 2 s e (fl,v,T,t)v 



WmW^t") 



(G3) 

(G4) 
(G5) 



a ,m.m' 
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Figure Fl. Probing the gaussianity of the harmonic expansions ai m , de- 
scribing the flux density of mass (Top) and bg m , describing the potential 
i/> e (bottom). Units are arbitrary. Only the real part of the coefficients is 
shown here, but the imaginary parts have similar distributions. Clearly, the 
distributions are quasi-gaussians. 



where 

K°v] m (t) = y^(3cr 2 /xq + fil)T m (t), 

a 

and 

/ 2ir 

T m (t) = J — {c^ tl _i{t) - c™ hl (t))e e 



(G6) 



— (-c"i _i(t) - c™ M (i))e,£ 



+2y gC^ li0 (t)e r . 



(G7) 



Using equation <G8t and equation <G9t . we find: 

w pVi (ft, t)zu pVj (CI, t) 



ttpOiOj (d, t) + 



U7 p (£l,t) 



di'dlY' s e (f2, v, r, t)viVj 



5^[q»(*)L^m(n). 



(G9) 



(G10) 



(Gil) 



The six independent elements of the symmetric tensor q(t) can be 
computed from c™ m , (t) coefficients using equation equation <G3t 
and recalling that: 



/ 



(2h + l)(2l 2 + l) (2l 3 + l) , ; 



*1 4 4 
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where 



li h 4 
mi rri2 mz 
symbol. One can find: 



mi 77i2 mz 
= WSifi&rn. is the Wigner 3-j 



fG12) 



(G13) 



am' 
am' 
am' 
am' 
am' 
am' 

where H™ m ,(t) = + l)(6a 2 ^ + 

APPENDIX H: NOTATIONS 



G2 From coefficients to advected velocity dispersion 

The distribution of advected velocity dispersion is given by : 

zu paitTj (n,v,r,t) = s e (n,v,r,t)(v-v(n,i)) i (v-v(n,t)) i ,(G8) 

where the subscripts i and j stand for r,9,cf> and where 
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Symbol Meaning 



(r, v) position and velocity 

t, t lookback time variables 

i?200 virial radius measured at z=0 

Vc circular velocity measured at -R200 

F the phase space distribution function of the halo 

ip the self-gravitating potential 

i/> e the external potential, induced by external perturbations 

s e the source function in phase space 

■uu x the flux density of x 

■uup the flux density of mass 

OTpv the flux density of momentum 

TOp CT(T the flux density of velocity dispersion 

i/il™] (r) 3D projection basis of the potential 

</>' n l (r) 6D projection basis of the source term 

a(t) expansion coefficients of the potential/density response 

b(t) expansion coefficients of the external potential perturbation 

c(t) expansion coefficients of the source perturbation 

f2 angular position on the virial sphere (2 angles) 

r angular orientation of the velocity vectore on the virial sphere (2 angles) 

v velocity's amplitude 

X angular average of X 

X temporal average of X 

(X) statistical expectation (or average value) of X 

(X) most probable value (or mode) of X 

5[x] (^) contrast density of X measured on the virial sphere 

m = (£, m) harmonic coefficients related to 

m' = (£' ,m') harmonic coefficients related to T 

a e.m(t) harmonic expansion coefficients of ] 

bl,m(t) harmonic expansion coefficients of Su,e] 

Ci (t) angular power spectrum 

T|(t, i + At) angular-temporal power spectrum 

w(6, t,t + At) angulo-temporal correlation function measured on the sphere for an angulo-temporal separation 6 and At 

<& M (t) accretion rate measured at the virial radius (averaged over all directions) 

■d(Ti,t) PDF of the velocity's incidence angle 

b impact parameter 

(p(v, t) PDF of the velocity's amplitude 

p(v, t) Joint PDF of the velocity's incidence angle and amplitude 



Table HI. A summary of the notations used throughout the paper. 



